added a new module called numerics.f90 which reads in all numerical "parameters" from the file numerics.config (also being added). From now on this file has to be located in the working directory of the FEM-model and has to contain all necessary parameters.

trunk/CPFEM.f90

@@ -1,221 +1,297 @@
 !##############################################################
- MODULE CPFEM
+MODULE CPFEM
 !##############################################################
 !    *** CPFEM engine ***
 !
- use prec, only: pReal,pInt
+use prec, only:                                   pReal, &
- implicit none
+                                                  pInt
-!
+implicit none
-! ****************************************************************
-! *** General variables for the material behaviour calculation ***
-! ****************************************************************
- real(pReal), dimension (:,:,:),        allocatable :: CPFEM_cs               ! Cauchy stress
- real(pReal), dimension (:,:,:,:),      allocatable :: CPFEM_dcsdE            ! Cauchy stress tangent
- real(pReal), dimension (:,:,:,:),      allocatable :: CPFEM_dcsdE_knownGood  ! known good tangent
  
- logical :: CPFEM_init_done    = .false.       ! remember whether init has been done already
+real(pReal), parameter ::                         CPFEM_odd_stress    = 1e15_pReal, &
- logical :: CPFEM_calc_done    = .false.       ! remember whether first IP has already calced the results
+                                                  CPFEM_odd_jacobian  = 1e50_pReal
+
+real(pReal), dimension (:,:,:),   allocatable ::  CPFEM_cs                            ! Cauchy stress
+real(pReal), dimension (:,:,:,:), allocatable ::  CPFEM_dcsdE                         ! Cauchy stress tangent
+real(pReal), dimension (:,:,:,:), allocatable ::  CPFEM_dcsdE_knownGood               ! known good tangent
+
+logical ::                                        CPFEM_init_done     = .false., &    ! remember whether init has been done already
+                                                  CPFEM_calc_done     = .false.       ! remember whether first IP has already calced the results
+ 
+
+CONTAINS
 
- real(pReal), parameter :: CPFEM_odd_stress = 1e15_pReal, CPFEM_odd_jacobian = 1e50_pReal
-!
- CONTAINS
-!
 !*********************************************************
 !***    allocate the arrays defined in module CPFEM    ***
 !***    and initialize them                            ***
 !*********************************************************
- SUBROUTINE CPFEM_init()
+subroutine CPFEM_init()
 
- use prec,           only: pInt,pReal
+  use prec, only:                                 pInt
- use FEsolving,      only: parallelExecution,symmetricSolver,FEsolving_execElem,FEsolving_execIP
+  use FEsolving, only:                            parallelExecution, &
- use mesh,           only: mesh_element,mesh_NcpElems,mesh_maxNips,FE_Nips
+                                                  symmetricSolver
- use material,       only: homogenization_maxNgrains
+  use mesh, only:                                 mesh_NcpElems, &
- use constitutive,   only: constitutive_maxSizePostResults
+                                                  mesh_maxNips
- use crystallite,    only: crystallite_Nresults
- use homogenization, only: homogenization_maxSizePostResults
 
- implicit none
+  implicit none
- integer(pInt) e,i,g
 
- allocate(CPFEM_cs(6,mesh_maxNips,mesh_NcpElems)) ;                CPFEM_cs              = 0.0_pReal
+  ! initialize stress and jacobian to zero 
- allocate(CPFEM_dcsdE(6,6,mesh_maxNips,mesh_NcpElems)) ;           CPFEM_dcsde           = 0.0_pReal
+  allocate(CPFEM_cs(6,mesh_maxNips,mesh_NcpElems)) ;                CPFEM_cs              = 0.0_pReal
- allocate(CPFEM_dcsdE_knownGood(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dcsde_knownGood = 0.0_pReal
+  allocate(CPFEM_dcsdE(6,6,mesh_maxNips,mesh_NcpElems)) ;           CPFEM_dcsde           = 0.0_pReal
+  allocate(CPFEM_dcsdE_knownGood(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dcsde_knownGood = 0.0_pReal
+
+  !$OMP CRITICAL (write2out)
+    write(6,*)
+    write(6,*) '<<<+-  cpfem init  -+>>>'
+    write(6,*)
+    write(6,'(a32,x,6(i5,x))') 'CPFEM_cs:              ', shape(CPFEM_cs)
+    write(6,'(a32,x,6(i5,x))') 'CPFEM_dcsde:           ', shape(CPFEM_dcsde)
+    write(6,'(a32,x,6(i5,x))') 'CPFEM_dcsde_knownGood: ', shape(CPFEM_dcsde_knownGood)
+    write(6,*)
+    write(6,*) 'parallelExecution:    ', parallelExecution
+    write(6,*) 'symmetricSolver:      ', symmetricSolver
+    call flush(6)
+  !$OMP END CRITICAL (write2out)
+  return
+
+endsubroutine
 
 
-!    *** Output to MARC output file ***
-!$OMP CRITICAL (write2out)
- write(6,*)
- write(6,*) '<<<+-  cpfem init  -+>>>'
- write(6,*)
- write(6,'(a32,x,6(i5,x))') 'CPFEM_cs:              ', shape(CPFEM_cs)
- write(6,'(a32,x,6(i5,x))') 'CPFEM_dcsde:           ', shape(CPFEM_dcsde)
- write(6,'(a32,x,6(i5,x))') 'CPFEM_dcsde_knownGood: ', shape(CPFEM_dcsde_knownGood)
- write(6,*)
- write(6,*) 'parallelExecution:    ', parallelExecution
- write(6,*) 'symmetricSolver:      ', symmetricSolver
- call flush(6)
-!$OMP END CRITICAL (write2out)
- return
-!
- END SUBROUTINE
-!
-!
 !***********************************************************************
 !***    perform initialization at first call, update variables and   ***
 !***    call the actual material model                               ***
-!
-!     CPFEM_mode             computation mode (regular, collection, recycle)
-!     ffn                    deformation gradient for t=t0
-!     ffn1                   deformation gradient for t=t1
-!     Temperature            temperature
-!     CPFEM_dt               time increment
-!     CPFEM_en               element number
-!     CPFEM_in               intergration point number
-!     CPFEM_stress           stress vector in Mandel notation
-!     CPFEM_updateJaco       flag to initiate computation of Jacobian
-!     CPFEM_jaco             jacobian in Mandel notation
-!     CPFEM_ngens            size of stress strain law
 !***********************************************************************
-subroutine CPFEM_general(CPFEM_mode, ffn, ffn1, Temperature, CPFEM_dt,&
+subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchyStress, jacobian, ngens)
-                         CPFEM_en, CPFEM_in, CPFEM_stress, CPFEM_updateJaco, CPFEM_jaco, CPFEM_ngens)
+  ! note: cauchyStress = Cauchy stress cs(6) and jacobian = Consistent tangent dcs/de
-! note: CPFEM_stress = Cauchy stress cs(6) and CPFEM_jaco = Consistent tangent dcs/de
-!
- use prec,         only: pReal,pInt
- use FEsolving
- use debug
- use math
- use mesh,         only: mesh_init,&
-                         mesh_FEasCP,mesh_element,mesh_NcpElems,mesh_maxNips,FE_Nips
- use lattice,      only: lattice_init
- use material,     only: material_init, homogenization_maxNgrains
- use constitutive, only: constitutive_init,&
-                         constitutive_state0,constitutive_state
- use crystallite
- use homogenization 
- implicit none
 
- integer(pInt) CPFEM_en, CPFEM_in, cp_en, CPFEM_ngens, i,j,k,l,m,n
+  !*** variables and functions from other modules ***!
- real(pReal), dimension (3,3)        :: ffn,ffn1,Kirchhoff
+  use prec, only:                                     pReal, &
- real(pReal), dimension (3,3,3,3)    :: H, H_sym
+                                                      pInt
- real(pReal), dimension(CPFEM_ngens) :: CPFEM_stress
+  use numerics, only:                                 numerics_init, & 
- real(pReal), dimension(CPFEM_ngens,CPFEM_ngens) :: CPFEM_jaco
+                                                      relevantStrain, &
- real(pReal) Temperature,CPFEM_dt,J_inverse
+                                                      iJacoStiffness
- integer(pInt) CPFEM_mode               ! 1: regular computation with aged results&
+  use debug, only:                                    debug_init
-                                        ! 2: regular computation&
+  use FEsolving, only:                                FE_init, &
-                                        ! 3: collection of FEM data&
+                                                      parallelExecution, &
-                                        ! 4: recycling of former results (MARC speciality)&
+                                                      outdatedFFN1, &
-                                        ! 5: record tangent from former converged inc&
+                                                      cycleCounter, &
-                                        ! 6: restore tangent from former converged inc
+                                                      theInc, &
- integer(pInt) e
+                                                      theCycle, &
- logical CPFEM_updateJaco
+                                                      theLovl, &
+                                                      theTime, &
+                                                      FEsolving_execElem, &
+                                                      FEsolving_execIP
+  use math, only:                                     math_init, &
+                                                      math_identity2nd, &
+                                                      math_mul33x33, &
+                                                      math_det3x3, &
+                                                      math_I3, &
+                                                      math_Mandel3333to66, &
+                                                      math_Mandel33to6
+  use mesh, only:                                     mesh_init, &
+                                                      mesh_FEasCP, &
+                                                      mesh_NcpElems, &
+                                                      mesh_maxNips
+  use lattice, only:                                  lattice_init
+  use material, only:                                 material_init, &
+                                                      homogenization_maxNgrains
+  use constitutive, only:                             constitutive_init,&
+                                                      constitutive_state0,constitutive_state
+  use crystallite, only:                              crystallite_init, &
+                                                      crystallite_F0, &
+                                                      crystallite_partionedF, &
+                                                      crystallite_Fp0, &
+                                                      crystallite_Fp, &
+                                                      crystallite_Lp0, &
+                                                      crystallite_Lp
+  use homogenization, only:                           homogenization_init, &
+                                                      homogenization_sizeState, &
+                                                      homogenization_state, &
+                                                      homogenization_state0, &
+                                                      materialpoint_F, &
+                                                      materialpoint_F0, &
+                                                      materialpoint_P, &
+                                                      materialpoint_dPdF, &
+                                                      materialpoint_Temperature, &
+                                                      materialpoint_stressAndItsTangent, &
+                                                      materialpoint_postResults
   
- if (.not. CPFEM_init_done) then        ! initialization step (three dimensional stress state check missing?)
+  implicit none
-   call math_init()
-   call FE_init()
-   call mesh_init()
   
-   FEsolving_execElem = (/1,mesh_NcpElems/)
+  !*** input variables ***!
-   allocate(FEsolving_execIP(2,mesh_NcpElems)); FEsolving_execIP = 1_pInt
+  integer(pInt), intent(in) ::                        element, &          ! FE element number
-   forall (e = 1:mesh_NcpElems) FEsolving_execIP(2,e) = FE_Nips(mesh_element(2,e))
+                                                      IP, &               ! FE integration point number
+                                                      ngens               ! size of stress strain law
+  real(pReal), intent(in) ::                          Temperature, &      ! temperature
+                                                      dt                  ! time increment
+  real(pReal), dimension (3,3), intent(in) ::         ffn, &              ! deformation gradient for t=t0
+                                                      ffn1                ! deformation gradient for t=t1
+  integer(pInt), intent(in) ::                        mode                ! computation mode  1: regular computation with aged results
+                                                                          !                   2: regular computation
+                                                                          !                   3: collection of FEM data
+                                                                          !                   4: recycling of former results (MARC speciality)
+                                                                          !                   5: record tangent from former converged inc
+                                                                          !                   6: restore tangent from former converged inc
   
-   call lattice_init()
+  !*** output variables ***!
-   call material_init()
+  real(pReal), dimension(ngens), intent(out) ::       cauchyStress        ! stress vector in Mandel notation
-   call constitutive_init()
+  real(pReal), dimension(ngens,ngens), intent(out) :: jacobian            ! jacobian in Mandel notation
-   call crystallite_init()
-   call homogenization_init()
-   call CPFEM_init()
-   CPFEM_init_done = .true.
- endif
   
- cp_en = mesh_FEasCP('elem',CPFEM_en)
+  !*** local variables ***!
- if (cp_en == 1 .and. CPFEM_in == 1) then
+  real(pReal)                                         J_inverse           ! inverse of Jacobian
+  real(pReal), dimension (3,3) ::                     Kirchhoff
+  real(pReal), dimension (3,3,3,3) ::                 H, &
+                                                      H_sym
+  integer(pInt)                                       cp_en, &            ! crystal plasticity element number
+                                                      i, &
+                                                      j, &
+                                                      k, &
+                                                      l, &
+                                                      m, &
+                                                      n
+  logical                                             updateJaco          ! flag indicating if JAcobian has to be updated
+  
+  !*** global variables ***!
+  ! CPFEM_cs, &
+  ! CPFEM_dcsdE, &
+  ! CPFEM_dcsdE_knownGood, &
+  ! CPFEM_init_done, &
+  ! CPFEM_calc_done, &
+  ! CPFEM_odd_stress, &
+  ! CPFEM_odd_jacobian
+  
+
+  ! initialization step (three dimensional stress state check missing?)
+  if (.not. CPFEM_init_done) then
+    call numerics_init()
+    call debug_init()
+    call math_init()
+    call FE_init()
+    call mesh_init()
+    call lattice_init()
+    call material_init()
+    call constitutive_init()
+    call crystallite_init()
+    call homogenization_init()
+    call CPFEM_init()
+    CPFEM_init_done = .true.
+  endif
+
+  cp_en = mesh_FEasCP('elem',element)
+  
+  if (cp_en == 1 .and. IP == 1) then
     write(6,*) '#####################################'
     write(6,'(a10,1x,f8.4,1x,a10,1x,i4,1x,a10,1x,i3,1x,a10,1x,i2,x,a10,1x,i2)') &
     'theTime',theTime,'theInc',theInc,'theCycle',theCycle,'theLovl',theLovl,&
-    'mode',CPFEM_mode
+    'mode',mode
     write(6,*) '#####################################'
- endif
+  endif
 
- select case (CPFEM_mode)
+  ! according to our "mode" we decide what to do
-    case (1,2)     ! regular computation (with aging of results if mode == 1)
+  select case (mode)
-       if (CPFEM_mode == 1) then                  ! age results at start of new increment
+    
-         crystallite_F0  = crystallite_partionedF ! crystallite deformation (_subF is perturbed...)
+    ! --+>> REGULAR COMPUTATION (WITH AGING OF RESULTS IF MODE == 1) <<+-- 
-         crystallite_Fp0 = crystallite_Fp         ! crystallite plastic deformation
+    case (1,2)
-         crystallite_Lp0 = crystallite_Lp         ! crystallite plastic velocity
+      ! age results if mode == 1
-         forall (i = 1:homogenization_maxNgrains,&
+      if (mode == 1) then
+        crystallite_F0  = crystallite_partionedF                          ! crystallite deformation (_subF is perturbed...)
+        crystallite_Fp0 = crystallite_Fp                                  ! crystallite plastic deformation
+        crystallite_Lp0 = crystallite_Lp                                  ! crystallite plastic velocity
+        forall ( i = 1:homogenization_maxNgrains, &
                  j = 1:mesh_maxNips, &
-                 k = 1:mesh_NcpElems) &
+                 k = 1:mesh_NcpElems ) &
-           constitutive_state0(i,j,k)%p = constitutive_state(i,j,k)%p   ! microstructure of crystallites
+          constitutive_state0(i,j,k)%p = constitutive_state(i,j,k)%p      ! microstructure of crystallites
-         write(6,'(a10,/,4(3(f10.3,x),/))') 'aged state',constitutive_state(1,1,1)%p/1e6
+        write(6,'(a10,/,4(3(f10.3,x),/))') 'aged state',constitutive_state(1,1,1)%p/1e6
-         do j = 1,mesh_maxNips
+        do j = 1,mesh_maxNips
-           do k = 1,mesh_NcpElems
+          do k = 1,mesh_NcpElems
-             if (homogenization_sizeState(j,k) > 0_pInt) &
+            if (homogenization_sizeState(j,k) > 0_pInt) &
-               homogenization_state0(j,k)%p = homogenization_state(j,k)%p   ! internal state of homogenization scheme
+              homogenization_state0(j,k)%p = homogenization_state(j,k)%p  ! internal state of homogenization scheme
-           enddo
+          enddo
-         enddo
+        enddo
-       endif
-
-       if (outdatedFFN1 .or. any(abs(ffn1 - materialpoint_F(:,:,CPFEM_in,cp_en)) > relevantStrain)) then
-         if (.not. outdatedFFN1) write(6,'(a11,x,i5,x,i2,x,a10,/,3(3(f10.3,x),/))') 'outdated at',cp_en,CPFEM_in,'FFN1 now:',ffn1(:,1),ffn1(:,2),ffn1(:,3)
-         outdatedFFN1 = .true.
-         CPFEM_cs(1:CPFEM_ngens,CPFEM_in,cp_en)                  = CPFEM_odd_stress
-         CPFEM_dcsde(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_jacobian*math_identity2nd(CPFEM_ngens)
-       else
-         if (.not. parallelExecution) then
-           FEsolving_execElem(1) = cp_en
-           FEsolving_execElem(2) = cp_en
-           FEsolving_execIP(1,cp_en) = CPFEM_in
-           FEsolving_execIP(2,cp_en) = CPFEM_in
-           call materialpoint_stressAndItsTangent(CPFEM_updateJaco, CPFEM_dt)
-           call materialpoint_postResults(CPFEM_dt)
-         elseif (.not. CPFEM_calc_done) then
-           call materialpoint_stressAndItsTangent(CPFEM_updateJaco, CPFEM_dt) ! parallel execution inside
-           call materialpoint_postResults(CPFEM_dt)
-           CPFEM_calc_done = .true.
-         endif
-
-!  translate from P and dP/dF to CS and dCS/dE
-         Kirchhoff = math_mul33x33(materialpoint_P(:,:,CPFEM_in, cp_en),transpose(materialpoint_F(:,:,CPFEM_in, cp_en)))
-         J_inverse  = 1.0_pReal/math_det3x3(materialpoint_F(:,:,CPFEM_in, cp_en))
-         CPFEM_cs(1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel33to6(J_inverse*Kirchhoff)
-
-         H = 0.0_pReal
-         forall(i=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3) &
-           H(i,j,k,l) = H(i,j,k,l) + &
-                          materialpoint_F(j,m,CPFEM_in,cp_en) * &
-                          materialpoint_F(l,n,CPFEM_in,cp_en) * &
-                          materialpoint_dPdF(i,m,k,n,CPFEM_in,cp_en) - &
-                          math_I3(j,l)*materialpoint_F(i,m,CPFEM_in,cp_en)*materialpoint_P(k,m,CPFEM_in,cp_en) + &
-                          0.5_pReal*(math_I3(i,k)*Kirchhoff(j,l) + math_I3(j,l)*Kirchhoff(i,k) + &
-                                     math_I3(i,l)*Kirchhoff(j,k) + math_I3(j,k)*Kirchhoff(i,l))
-         forall(i=1:3,j=1:3,k=1:3,l=1:3) &
-            H_sym(i,j,k,l)= 0.25_pReal*(H(i,j,k,l)+H(j,i,k,l)+H(i,j,l,k)+H(j,i,l,k))          ! where to use the symmetric version??
-         CPFEM_dcsde(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel3333to66(J_inverse*H)
       endif
-    case (3)    ! collect and return odd result
-       materialpoint_Temperature(CPFEM_in,cp_en) = Temperature
-       materialpoint_F0(:,:,CPFEM_in,cp_en)      = ffn
-       materialpoint_F(:,:,CPFEM_in,cp_en)       = ffn1
-       CPFEM_cs(1:CPFEM_ngens,CPFEM_in,cp_en)                  = CPFEM_odd_stress
-       CPFEM_dcsde(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_jacobian*math_identity2nd(CPFEM_ngens)
-       CPFEM_calc_done = .false.
 
-    case (4)    ! do nothing since we can recycle the former results (MARC specialty)
+      ! deformation gradient outdated or any actual deformation gradient differs more than relevantStrain from the stored one
-    case (5)    ! record consistent tangent at beginning of new FE increment (while recycling)
+      if (outdatedFFN1 .or. any(abs(ffn1 - materialpoint_F(:,:,IP,cp_en)) > relevantStrain)) then
+        if (.not. outdatedFFN1) & 
+          write(6,'(a11,x,i5,x,i2,x,a10,/,3(3(f10.3,x),/))') 'outdated at',cp_en,IP,'FFN1 now:',ffn1(:,1),ffn1(:,2),ffn1(:,3)
+        outdatedFFN1 = .true.
+        CPFEM_cs(1:ngens,IP,cp_en)              = CPFEM_odd_stress
+        CPFEM_dcsde(1:ngens,1:ngens,IP,cp_en)   = CPFEM_odd_jacobian*math_identity2nd(ngens)
+      
+      ! deformation gradient is not outdated
+      else
+        ! set flag for Jacobian update
+        updateJaco = (mod(cycleCounter-4,4_pInt*iJacoStiffness)==0)
+        
+        ! no parallel computation
+        if (.not. parallelExecution) then
+          ! we just take one single element and IP
+          FEsolving_execElem(1)     = cp_en
+          FEsolving_execElem(2)     = cp_en
+          FEsolving_execIP(1,cp_en) = IP
+          FEsolving_execIP(2,cp_en) = IP
+          call materialpoint_stressAndItsTangent(updateJaco, dt)          ! calculate stress and its tangent
+          call materialpoint_postResults(dt)                              ! post results
+          
+        ! parallel computation and calulation not yet done
+        elseif (.not. CPFEM_calc_done) then
+          call materialpoint_stressAndItsTangent(updateJaco, dt)          ! calculate stress and its tangent (parallel execution inside)
+          call materialpoint_postResults(dt)                              ! post results
+          CPFEM_calc_done = .true.
+        endif
+        
+        !  translate from P to CS
+        Kirchhoff = math_mul33x33(materialpoint_P(:,:,IP, cp_en),transpose(materialpoint_F(:,:,IP, cp_en)))
+        J_inverse  = 1.0_pReal/math_det3x3(materialpoint_F(:,:,IP, cp_en))
+        CPFEM_cs(1:ngens,IP,cp_en) = math_Mandel33to6(J_inverse*Kirchhoff)
+
+        !  translate from dP/dF to dCS/dE
+        H = 0.0_pReal
+        forall(i=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3) &
+          H(i,j,k,l) = H(i,j,k,l) + &
+                        materialpoint_F(j,m,IP,cp_en) * &
+                        materialpoint_F(l,n,IP,cp_en) * &
+                        materialpoint_dPdF(i,m,k,n,IP,cp_en) - &
+                        math_I3(j,l)*materialpoint_F(i,m,IP,cp_en)*materialpoint_P(k,m,IP,cp_en) + &
+                        0.5_pReal*(math_I3(i,k)*Kirchhoff(j,l) + math_I3(j,l)*Kirchhoff(i,k) + &
+                                   math_I3(i,l)*Kirchhoff(j,k) + math_I3(j,k)*Kirchhoff(i,l))
+        forall(i=1:3,j=1:3,k=1:3,l=1:3) &
+          H_sym(i,j,k,l)= 0.25_pReal*(H(i,j,k,l)+H(j,i,k,l)+H(i,j,l,k)+H(j,i,l,k))  ! where to use the symmetric version??
+        CPFEM_dcsde(1:ngens,1:ngens,IP,cp_en) = math_Mandel3333to66(J_inverse*H)
+        
+      endif
+    
+    ! --+>> COLLECTION OF FEM DATA AND RETURN OF ODD STRESS AND JACOBIAN <<+-- 
+    case (3)
+      materialpoint_Temperature(IP,cp_en)   = Temperature
+      materialpoint_F0(:,:,IP,cp_en)        = ffn
+      materialpoint_F(:,:,IP,cp_en)         = ffn1
+      CPFEM_cs(1:ngens,IP,cp_en)            = CPFEM_odd_stress
+      CPFEM_dcsde(1:ngens,1:ngens,IP,cp_en) = CPFEM_odd_jacobian*math_identity2nd(ngens)
+      CPFEM_calc_done = .false.
+    
+    ! --+>> RECYCLING OF FORMER RESULTS (MARC SPECIALTY) <<+--
+    case (4)
+      ! do nothing
+    
+    ! --+>> RECORD JACOBIAN FROM FORMER CONVERGED INC <<+--
+    case (5)
        CPFEM_dcsde_knownGood = CPFEM_dcsde
-    case (6)    ! restore consistent tangent after FE cutback
+    
+    ! --+>> RESTORE CONSISTENT JACOBIAN FROM FORMER CONVERGED INC <<+--
+    case (6)
        CPFEM_dcsde = CPFEM_dcsde_knownGood
- end select
        
-! return the local stress and the jacobian from storage
+  end select
- CPFEM_stress(1:CPFEM_ngens) = CPFEM_cs(1:CPFEM_ngens,CPFEM_in,cp_en)
- CPFEM_jaco(1:CPFEM_ngens,1:CPFEM_ngens) = CPFEM_dcsdE(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en)
 
- return
+  ! return the local stress and the jacobian from storage
+  cauchyStress(1:ngens) = CPFEM_cs(1:ngens,IP,cp_en)
+  jacobian(1:ngens,1:ngens) = CPFEM_dcsdE(1:ngens,1:ngens,IP,cp_en)
+
+  return
 
 end subroutine
 
 
 
- END MODULE
+ END MODULE CPFEM
-!##############################################################
+ 

trunk/IO.f90

@@ -30,7 +30,7 @@
 ! open existing file to given unit
 ! path to file is relative to working directory
 !********************************************************************
- logical FUNCTION IO_open_file(unit,relPath)
+ logical function IO_open_file(unit,relPath)
 
  use prec, only: pInt
  implicit none
@@ -47,13 +47,13 @@
 100 IO_open_file = .false.
  return
  
- END FUNCTION
+ endfunction
 
 
 !********************************************************************
 ! open FEM inputfile to given unit
 !********************************************************************
- logical FUNCTION IO_open_inputFile(unit)
+ logical function IO_open_inputFile(unit)
 
  use prec, only: pReal, pInt
  implicit none
@@ -68,20 +68,20 @@
      ext='dat' ! MARC
  else
      ext='inp' ! ABAQUS
- end if
+ endif
  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
+ endfunction
 
 
 !********************************************************************
 ! hybrid IA repetition counter
 !********************************************************************
- FUNCTION hybridIA_reps(dV_V,steps,C)
+ function hybridIA_reps(dV_V,steps,C)
 
  use prec, only: pReal, pInt
  implicit none
@@ -96,18 +96,18 @@
    do Phi =1,steps(2)
      do phi2=1,steps(3)
        hybridIA_reps = hybridIA_reps+nint(C*dV_V(phi2,Phi,phi1), pInt)
-     end do
+     enddo
-   end do
+   enddo
- end do
+ enddo
  return
  
- END FUNCTION
+ endfunction
 
 
 !********************************************************************
 ! hybrid IA sampling of ODFfile
 !********************************************************************
- FUNCTION IO_hybridIA(Nast,ODFfileName)
+ function IO_hybridIA(Nast,ODFfileName)
 
  use prec, only: pReal, pInt
  use math, only: inRad
@@ -135,7 +135,7 @@
  if (pos(1).ne.3) goto 100
  do i=1,3
    limits(i) = IO_intValue(line,pos,i)*inRad
- end do
+ enddo
 
 !--- deltas in phi1, Phi, phi2 ---
  read(999,fmt=fileFormat,end=100) line
@@ -143,7 +143,7 @@
  if (pos(1).ne.3) goto 100
  do i=1,3
    deltas(i) = IO_intValue(line,pos,i)*inRad
- end do
+ enddo
  steps = nint(limits/deltas,pInt)
  allocate(dV_V(steps(3),steps(2),steps(1)))
 
@@ -153,7 +153,7 @@
    center = 0.5_pReal
  else
    center = 0.0_pReal
- end if
+ endif
  
 !--- skip blank line ---
  read(999,fmt=fileFormat,end=100) line
@@ -172,11 +172,11 @@
          sum_dV_V = sum_dV_V+prob
        else
          prob = 0.0_pReal
-       end if
+       endif
        dV_V(phi2,Phi,phi1) = prob*dg_0*sin((Phi-1.0_pReal+center)*deltas(2))
-     end do
+     enddo
-   end do
+   enddo
- end do  
+ enddo  
  
  dV_V = dV_V/sum_dV_V  ! normalize to 1
  
@@ -188,7 +188,7 @@
  do while (hybridIA_reps(dV_V,steps,upperC) < Nset)
    lowerC = upperC
    upperC = upperC*2.0_pReal
- end do
+ enddo
 !--- binary search for best C ---
  do
    C = (upperC+lowerC)/2.0_pReal
@@ -203,8 +203,8 @@
      upperC = C
    else
      exit
-   end if
+   endif
- end do
+ enddo
  
  allocate(binSet(Nreps))
  bin = 0 ! bin counter
@@ -216,9 +216,9 @@
        binSet(i:i+reps-1) = bin
        bin = bin+1 ! advance bin
        i = i+reps ! advance set
-     end do
+     enddo
-   end do
+   enddo
- end do
+ enddo
 
  do i=1,Nast
    if (i < Nast) then
@@ -226,13 +226,13 @@
      j = nint(rnd*(Nast-i)+i+0.5_pReal,pInt)
    else
      j = i
-   end if
+   endif
    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
+ enddo
  close(999)
  return
 
@@ -241,13 +241,13 @@
  close(999)
  return
  
- END FUNCTION 
+ endfunction 
 
 
 !********************************************************************
 ! identifies lines without content
 !********************************************************************
- PURE FUNCTION IO_isBlank (line)
+ pure function IO_isBlank (line)
 
  use prec, only: pInt
  implicit none
@@ -264,12 +264,12 @@
 
  return
  
- END FUNCTION
+ endfunction
 
 !********************************************************************
 ! get tagged content of line
 !********************************************************************
- PURE FUNCTION IO_getTag (line,openChar,closechar)
+ pure function IO_getTag (line,openChar,closechar)
 
  use prec, only: pInt
  implicit none
@@ -288,11 +288,11 @@
 
  return
  
- END FUNCTION
+ endfunction
 
 
 !*********************************************************************
- FUNCTION IO_countSections(file,part)
+ function IO_countSections(file,part)
 !*********************************************************************
  use prec, only: pInt
  implicit none
@@ -321,13 +321,13 @@
 
 100 return
 
- END FUNCTION
+ endfunction
  
 
 !*********************************************************************
 ! return array of myTag counts within <part> for at most N[sections]
 !*********************************************************************
- FUNCTION IO_countTagInPart(file,part,myTag,Nsections)
+ function IO_countTagInPart(file,part,myTag,Nsections)
 
  use prec, only: pInt
  implicit none
@@ -367,13 +367,13 @@
 100 IO_countTagInPart = counter
  return
 
-END FUNCTION
+endfunction
 
 
 !*********************************************************************
 ! return array of myTag presence within <part> for at most N[sections]
 !*********************************************************************
- FUNCTION IO_spotTagInPart(file,part,myTag,Nsections)
+ function IO_spotTagInPart(file,part,myTag,Nsections)
 
  use prec, only: pInt
  implicit none
@@ -412,7 +412,7 @@ END FUNCTION
 
 100 return
 
-END FUNCTION
+endfunction
 
 
 !********************************************************************
@@ -420,7 +420,7 @@ END FUNCTION
 ! return array containing number of parts found and
 ! their left/right positions to be used by IO_xxxVal
 !********************************************************************
- PURE FUNCTION IO_stringPos (line,N)
+ pure function IO_stringPos (line,N)
 
  use prec, only: pReal,pInt
  implicit none
@@ -438,17 +438,17 @@ END FUNCTION
    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
+ enddo
  IO_stringPos(1) = part-1
  return
 
- END FUNCTION
+ endfunction
 
 
 !********************************************************************
 ! read string value at pos from line
 !********************************************************************
- PURE FUNCTION IO_stringValue (line,positions,pos)
+ pure function IO_stringValue (line,positions,pos)
  
  use prec, only: pReal,pInt
  implicit none
@@ -464,13 +464,13 @@ END FUNCTION
  endif
  return
 
- END FUNCTION
+ endfunction
 
 
 !********************************************************************
 ! read string value at pos from fixed format line
 !********************************************************************
- PURE FUNCTION IO_fixedStringValue (line,ends,pos)
+ pure function IO_fixedStringValue (line,ends,pos)
  
  use prec, only: pReal,pInt
  implicit none
@@ -482,13 +482,13 @@ END FUNCTION
  IO_fixedStringValue = line(ends(pos)+1:ends(pos+1))
  return
 
- END FUNCTION
+ endfunction
 
 
 !********************************************************************
 ! read float value at pos from line
 !********************************************************************
- PURE FUNCTION IO_floatValue (line,positions,pos)
+ pure function IO_floatValue (line,positions,pos)
  
  use prec, only: pReal,pInt
  implicit none
@@ -504,13 +504,13 @@ END FUNCTION
 100 IO_floatValue = huge(1.0_pReal)
  return
 
- END FUNCTION
+ endfunction
 
 
 !********************************************************************
 ! read float value at pos from fixed format line
 !********************************************************************
- PURE FUNCTION IO_fixedFloatValue (line,ends,pos)
+ pure function IO_fixedFloatValue (line,ends,pos)
  
  use prec, only: pReal,pInt
  implicit none
@@ -524,13 +524,13 @@ END FUNCTION
 100 IO_fixedFloatValue = huge(1.0_pReal)
  return
 
- END FUNCTION
+ endfunction
 
 
 !********************************************************************
 ! read float x.y+z value at pos from format line line
 !********************************************************************
- PURE FUNCTION IO_fixedNoEFloatValue (line,ends,pos)
+ pure function IO_fixedNoEFloatValue (line,ends,pos)
  
  use prec, only: pReal,pInt
  implicit none
@@ -553,13 +553,13 @@ END FUNCTION
 100 IO_fixedNoEFloatValue = huge(1.0_pReal)
  return
 
- END FUNCTION
+ endfunction
 
 
 !********************************************************************
 ! read int value at pos from line
 !********************************************************************
- PURE FUNCTION IO_intValue (line,positions,pos)
+ pure function IO_intValue (line,positions,pos)
  
  use prec, only: pReal,pInt
  implicit none
@@ -575,13 +575,13 @@ END FUNCTION
 100 IO_intValue = huge(1_pInt)
  return
 
- END FUNCTION
+ endfunction
 
 
 !********************************************************************
 ! read int value at pos from fixed format line
 !********************************************************************
- PURE FUNCTION IO_fixedIntValue (line,ends,pos)
+ pure function IO_fixedIntValue (line,ends,pos)
  
  use prec, only: pReal,pInt
  implicit none
@@ -595,13 +595,13 @@ END FUNCTION
 100 IO_fixedIntValue = huge(1_pInt)
  return
 
- END FUNCTION
+ endfunction
 
 
 !********************************************************************
 ! change character in line to lower case
 !********************************************************************
- PURE FUNCTION IO_lc (line)
+ pure function IO_lc (line)
 
  use prec, only: pInt
  implicit none
@@ -616,13 +616,13 @@ END FUNCTION
  enddo
  return 
 
- END FUNCTION
+ endfunction
 
 
 !********************************************************************
 ! in place change of character in line to lower case
 !********************************************************************
- SUBROUTINE IO_lcInplace (line)
+ subroutine IO_lcInplace (line)
 
  use prec, only: pInt
  implicit none
@@ -638,13 +638,13 @@ END FUNCTION
  line = IO_lc
  return 
 
- END SUBROUTINE
+ endsubroutine
 
 
 !********************************************************************
 ! read on in file to skip (at least) N chunks (may be over multiple lines)
 !********************************************************************
- SUBROUTINE IO_skipChunks (unit,N)
+ subroutine IO_skipChunks (unit,N)
 
  use prec, only: pReal,pInt
  implicit none
@@ -659,17 +659,17 @@ END FUNCTION
    read(unit,'(A300)',end=100) line
    pos = IO_stringPos(line,maxNchunks)
    remainingChunks = remainingChunks - pos(1)
- end do
+ enddo
 100 return
 
- END SUBROUTINE
+ endsubroutine
 
  
 !********************************************************************
 ! count items in consecutive lines of ints concatenated by "c"
 ! as last char or range of values a "to" b
 !********************************************************************
- FUNCTION IO_countContinousIntValues (unit)
+ function IO_countContinousIntValues (unit)
 
  use prec, only: pReal,pInt
  implicit none
@@ -695,13 +695,13 @@ END FUNCTION
  enddo
 100 return
 
- END FUNCTION
+ endfunction
 
 !*********************************************************************
 ! read consecutive lines of ints concatenated by "c" as last char
 ! or range of values a "to" b
 !*********************************************************************
- FUNCTION IO_continousIntValues (unit,maxN,lookupName,lookupMap,lookupMaxN)
+ function IO_continousIntValues (unit,maxN,lookupName,lookupMap,lookupMaxN)
 
  use prec, only: pReal,pInt
  implicit none
@@ -746,7 +746,7 @@ END FUNCTION
  enddo
 100 return
 
- END FUNCTION
+ endfunction
 
 
 !********************************************************************
@@ -754,11 +754,9 @@ END FUNCTION
 ! and terminate the Marc run with exit #9xxx
 ! in ABAQUS either time step is reduced or execution terminated
 !********************************************************************
- SUBROUTINE IO_error(ID,e,i,g,ext_msg)
+ subroutine IO_error(ID,e,i,g,ext_msg)
 
  use prec, only: pInt
-
- use debug
  implicit none
 
  integer(pInt), intent(in) :: ID
@@ -819,6 +817,38 @@ END FUNCTION
    msg = 'Negative diffusion constant'
  case (240)
    msg = 'Non-positive Taylor factor'
+ case (260)
+   msg = 'Non-positive relevant strain'
+ case (261)
+   msg = 'Frequency for Stiffness update smaller than zero'
+ case (262)
+   msg = 'Frequency for Jacobian update of Lp residuum smaller than zero'
+ case (263)
+   msg = 'Non-positive perturbation value'
+ case (264)
+   msg = 'Limit for homogenization loop too small'
+ case (265)
+   msg = 'Limit for crystallite loop too small'
+ case (266)
+   msg = 'Limit for state loop too small'
+ case (267)
+   msg = 'Limit for stress loop too small'
+ case (268)
+   msg = 'Non-positive minimum substep size'
+ case (269)
+   msg = 'Non-positive relative tolerance for state'
+ case (270)
+   msg = 'Non-positive relative tolerance for stress'
+ case (271)
+   msg = 'Non-positive absolute tolerance for stress'
+ case (272)
+   msg = 'Non-positive relative tolerance of residual in GIA iteration'
+ case (273)
+   msg = 'Non-positive absolute tolerance of residual in GIA iteration'
+ case (274)
+   msg = 'Non-positive relative maximum value (upper bound) for GIA residual'
+ case (275)
+   msg = 'Limit for GIA iteration too small' 
  case (300)
    msg = 'This material can only be used with elements with three direct stress components'
  case (500)
@@ -851,7 +881,6 @@ END FUNCTION
  endif
  write(6,'(a38)') '+------------------------------------+'
  
- call debug_info()
  call flush(6)
  call quit(9000+ID)
 !$OMP END CRITICAL (write2out)
@@ -859,16 +888,15 @@ END FUNCTION
 ! ABAQUS returns in some cases
  return
 
- END SUBROUTINE
+ endsubroutine
 
 
 !********************************************************************
 ! write warning statements to standard out
 !********************************************************************
- SUBROUTINE IO_warning(ID,e,i,g,ext_msg)
+ subroutine IO_warning(ID,e,i,g,ext_msg)
 
  use prec, only: pInt
- use debug
  implicit none
 
  integer(pInt), intent(in) :: ID
@@ -901,6 +929,6 @@ END FUNCTION
  endif
  write(6,'(a38)') '+------------------------------------+'
 
- END SUBROUTINE
+ endsubroutine
  
  END MODULE IO

trunk/constitutive_j2.f90

@@ -176,7 +176,7 @@ subroutine constitutive_j2_init(file)
 
  return
 
-end subroutine
+endsubroutine
 
 
 function constitutive_j2_stateInit(myInstance)
@@ -193,7 +193,7 @@ function constitutive_j2_stateInit(myInstance)
  constitutive_j2_stateInit = constitutive_j2_s0(myInstance)
 
  return
-end function
+endfunction
 
 
 function constitutive_j2_homogenizedC(state,ipc,ip,el)
@@ -221,7 +221,7 @@ function constitutive_j2_homogenizedC(state,ipc,ip,el)
 
  return
 
-end function
+endfunction
 
 
 subroutine constitutive_j2_microstructure(Temperature,state,ipc,ip,el)
@@ -245,7 +245,7 @@ subroutine constitutive_j2_microstructure(Temperature,state,ipc,ip,el)
 
  matID = phase_constitutionInstance(material_phase(ipc,ip,el))
 
-end subroutine
+endsubroutine
 
 
 subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,state,ipc,ip,el)
@@ -310,7 +310,7 @@ subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,sta
  
  return
 
-end subroutine
+endsubroutine
 
 
 function constitutive_j2_dotState(Tstar_v,Temperature,state,ipc,ip,el)
@@ -351,7 +351,7 @@ function constitutive_j2_dotState(Tstar_v,Temperature,state,ipc,ip,el)
 
  return
 
-end function
+endfunction
 
 
 pure function constitutive_j2_postResults(Tstar_v,Temperature,dt,state,ipc,ip,el)
@@ -401,6 +401,6 @@ pure function constitutive_j2_postResults(Tstar_v,Temperature,dt,state,ipc,ip,el
  
  return
 
-end function
+endfunction
 
 END MODULE

trunk/constitutive_phenomenological.f90

@@ -180,7 +180,6 @@ subroutine constitutive_phenomenological_init(file)
 100 do i = 1,maxNinstance
 
    constitutive_phenomenological_structure(i) = lattice_initializeStructure(constitutive_phenomenological_structureName(i), &
-
                                                                       constitutive_phenomenological_CoverA(i))                                        ! sanity checks
    if (constitutive_phenomenological_structure(i) < 1 .or. &
        constitutive_phenomenological_structure(i) > 3)           call IO_error(201)

trunk/crystallite.f90

@@ -229,8 +229,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
 
  !*** variables and functions from other modules ***!
  use prec, only:                      pInt, &
-                                      pReal, &
+                                      pReal
-                                      subStepMin, &
+ use numerics, only:                  subStepMin, &
                                       pert_Fg, &
                                       nState, &
                                       nCryst
@@ -347,6 +347,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
  enddo
  !$OMPEND PARALLEL DO
  
+
  ! --+>> crystallite loop <<+--
 
  NiterationCrystallite = 0_pInt
@@ -458,7 +459,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
      !          constitutive_state is internally interpolated with .._subState0
      !          to account for substepping within _integrateStress
      ! results in crystallite_Fp,.._Lp
-     
      !$OMP PARALLEL DO
      do e = FEsolving_execElem(1),FEsolving_execElem(2)             ! iterate over elements to be processed
        myNgrains = homogenization_Ngrains(mesh_element(3,e))
@@ -610,7 +610,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
    enddo
    !$OMPEND PARALLEL DO
 
-  if (debugger) write (6,*) 'Stiffness calculation finished'
+   if (debugger) write (6,*) 'Stiffness calculation finished'
   
  endif
  
@@ -631,8 +631,8 @@ endsubroutine
  !*** variables and functions from other modules ***!
  use prec, only:                      pReal, &
                                       pInt, &
-                                      pLongInt, &
+                                      pLongInt
-                                      rTol_crystalliteState
+ use numerics, only:                  rTol_crystalliteState
  use constitutive, only:              constitutive_dotState, &
                                       constitutive_sizeDotState, &
                                       constitutive_subState0, &
@@ -708,12 +708,12 @@ endsubroutine
  !*** variables and functions from other modules ***!
  use prec, only:                      pReal, &
                                       pInt, &
-                                      pLongInt, &
+                                      pLongInt
-                                      nStress, &
+ use numerics, only:                  nStress, &
                                       aTol_crystalliteStress, &
                                       rTol_crystalliteStress, &
-                                      relevantStrain, &
+                                      iJacoLpresiduum, &
-                                      iJacoLpresiduum
+                                      relevantStrain
  use debug, only:                     debugger, &
                                       debug_cumLpCalls, &
                                       debug_cumLpTicks, &

trunk/debug.f90

@@ -4,12 +4,11 @@
 !##############################################################
  use prec
 
-
  implicit none
- integer(pInt), dimension(nStress) :: debug_StressLoopDistribution = 0_pInt
+ integer(pInt), dimension(:), allocatable :: debug_StressLoopDistribution
- integer(pInt), dimension(nState) :: debug_StateLoopDistribution = 0_pInt
+ integer(pInt), dimension(:), allocatable :: debug_StateLoopDistribution
- integer(pInt), dimension(nState) :: debug_StiffnessStateLoopDistribution = 0_pInt
+ integer(pInt), dimension(:), allocatable :: debug_StiffnessStateLoopDistribution
- integer(pInt), dimension(nCryst) :: debug_CrystalliteLoopDistribution = 0_pInt
+ integer(pInt), dimension(:), allocatable :: debug_CrystalliteLoopDistribution
  integer(pLongInt) :: debug_cumLpTicks = 0_pInt
  integer(pLongInt) :: debug_cumDotStateTicks = 0_pInt
  integer(pInt) :: debug_cumLpCalls = 0_pInt
@@ -19,36 +18,54 @@
 
  CONTAINS
 
+subroutine debug_init()
+  
+  use prec,     only: pInt  
+  use numerics, only: nStress, &
+                      nState, &
+                      nCryst
+  implicit none
+ 
+  allocate(debug_StressLoopDistribution(nStress)) ;        debug_StressLoopDistribution = 0_pInt
+  allocate(debug_StateLoopDistribution(nState)) ;          debug_StateLoopDistribution = 0_pInt
+  allocate(debug_StiffnessStateLoopDistribution(nState)) ; debug_StiffnessStateLoopDistribution = 0_pInt
+  allocate(debug_CrystalliteLoopDistribution(nCryst)) ;    debug_CrystalliteLoopDistribution = 0_pInt
+endsubroutine
+ 
 !********************************************************************
 ! reset debug distributions
 !********************************************************************
-SUBROUTINE debug_reset()
+subroutine debug_reset()
 
- use prec
+  use prec
- implicit none
+  implicit none
 
- debug_StressLoopDistribution         = 0_pInt ! initialize debugging data
+  debug_StressLoopDistribution         = 0_pInt ! initialize debugging data
- debug_StateLoopDistribution          = 0_pInt
+  debug_StateLoopDistribution          = 0_pInt
- debug_StiffnessStateLoopDistribution = 0_pInt
+  debug_StiffnessStateLoopDistribution = 0_pInt
- debug_CrystalliteLoopDistribution    = 0_pInt
+  debug_CrystalliteLoopDistribution    = 0_pInt
- debug_cumLpTicks       = 0_pInt
+  debug_cumLpTicks       = 0_pInt
- debug_cumDotStateTicks = 0_pInt
+  debug_cumDotStateTicks = 0_pInt
- debug_cumLpCalls       = 0_pInt
+  debug_cumLpCalls       = 0_pInt
- debug_cumDotStateCalls = 0_pInt
+  debug_cumDotStateCalls = 0_pInt
 
-END SUBROUTINE
+endsubroutine
 
 !********************************************************************
 ! write debug statements to standard out
 !********************************************************************
- SUBROUTINE debug_info()
+ subroutine debug_info()
 
  use prec
+ use numerics, only: nStress, &
+                      nState, &
+                      nCryst
  implicit none
 
  integer(pInt) i,integral
  integer(pLongInt) tickrate
  
+
  write(6,*)
  write(6,*) 'DEBUG Info'
  write(6,*)
@@ -112,6 +129,6 @@ END SUBROUTINE
  write(6,'(a15,i10,i10)') '          total',sum(debug_CrystalliteLoopDistribution),integral
  write(6,*)
 
- END SUBROUTINE
+ endsubroutine
  
  END MODULE debug

trunk/homogenization.f90

@@ -155,7 +155,7 @@ subroutine homogenization_init()
 
  return
 
-end subroutine
+endsubroutine
 
 
 !********************************************************************
@@ -168,11 +168,17 @@ subroutine materialpoint_stressAndItsTangent(&
      dt &             ! time increment
     )
 
- use prec,         only: pInt,pReal, subStepMin,nHomog
+ use prec, only:          pInt, &
- use FEsolving,    only: FEsolving_execElem, FEsolving_execIP
+                          pReal
- use mesh,         only: mesh_element
+ use numerics, only:      subStepMin, &
- use material,     only: homogenization_Ngrains
+                          nHomog
- use constitutive, only: constitutive_state0, constitutive_partionedState0, constitutive_state
+ use FEsolving, only:     FEsolving_execElem, &
+                          FEsolving_execIP
+ use mesh, only:          mesh_element
+ use material, only:      homogenization_Ngrains
+ use constitutive, only:  constitutive_state0, &
+                          constitutive_partionedState0, &
+                          constitutive_state
  use crystallite
  implicit none
  
@@ -331,7 +337,7 @@ subroutine materialpoint_stressAndItsTangent(&
  ! how to deal with stiffness?
  return
  
-end subroutine
+endsubroutine
 
 
 !********************************************************************
@@ -371,7 +377,7 @@ subroutine materialpoint_postResults(dt)
    enddo
 !$OMP END PARALLEL DO
 
- end subroutine
+ endsubroutine
  
  
 !********************************************************************
@@ -400,7 +406,7 @@ subroutine homogenization_partitionDeformation(&
                                                         homogenization_state(ip,el),ip,el)
  end select
 
-end subroutine
+endsubroutine
 
 
 !********************************************************************
@@ -431,7 +437,7 @@ function homogenization_updateState(&
 
  return
  
-end function
+endfunction
 
 
 !********************************************************************
@@ -459,7 +465,7 @@ subroutine homogenization_averageStressAndItsTangent(&
 
  return
  
-end subroutine
+endsubroutine
 
 
 !********************************************************************
@@ -489,6 +495,6 @@ function homogenization_postResults(&
 
  return
 
-end function
+endfunction
 
 END MODULE

trunk/homogenization_isostrain.f90

@@ -109,7 +109,7 @@ subroutine homogenization_isostrain_init(&
 
  return
 
-end subroutine
+endsubroutine
 
 
 !*********************************************************************
@@ -127,7 +127,7 @@ function homogenization_isostrain_stateInit(myInstance)
 
  return
  
-end function
+endfunction
 
 
 !********************************************************************
@@ -161,7 +161,7 @@ subroutine homogenization_isostrain_partitionDeformation(&
 
  return
 
-end subroutine
+endsubroutine
 
 
 !********************************************************************
@@ -195,7 +195,7 @@ function homogenization_isostrain_updateState(&
 
  return 
  
-end function
+endfunction
 
 
 !********************************************************************
@@ -232,7 +232,7 @@ subroutine homogenization_isostrain_averageStressAndItsTangent(&
 
  return
 
-end subroutine
+endsubroutine
 
 
 !********************************************************************
@@ -270,6 +270,6 @@ pure function homogenization_isostrain_postResults(&
  
  return
 
-end function
+endfunction
 
 END MODULE

trunk/lattice.f90

@@ -585,38 +585,38 @@ subroutine lattice_init()
 !**************************************
 !*      Module initialization         *
 !**************************************
- use IO, only: IO_open_file,IO_countSections,IO_countTagInPart,IO_error
+  use IO, only: IO_open_file,IO_countSections,IO_countTagInPart,IO_error
- use material, only: material_configfile,material_partPhase
+  use material, only: material_configfile,material_partPhase
- implicit none
+  implicit none
  
- integer(pInt), parameter :: fileunit = 200
+  integer(pInt), parameter :: fileunit = 200
- integer(pInt) i,Nsections
+  integer(pInt) i,Nsections
 
- if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file
+  if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file
- Nsections = IO_countSections(fileunit,material_partPhase)
+  Nsections = IO_countSections(fileunit,material_partPhase)
- lattice_Nstructure = 2_pInt + sum(IO_countTagInPart(fileunit,material_partPhase,'covera_ratio',Nsections)) ! fcc + bcc + all hex
+  lattice_Nstructure = 2_pInt + sum(IO_countTagInPart(fileunit,material_partPhase,'covera_ratio',Nsections)) ! fcc + bcc + all hex
- close(fileunit)
+  close(fileunit)
 
- allocate(lattice_Sslip(3,3,lattice_maxNslip,lattice_Nstructure)); lattice_Sslip = 0.0_pReal
+  allocate(lattice_Sslip(3,3,lattice_maxNslip,lattice_Nstructure)); lattice_Sslip = 0.0_pReal
- allocate(lattice_Sslip_v(6,lattice_maxNslip,lattice_Nstructure)); lattice_Sslip_v = 0.0_pReal
+  allocate(lattice_Sslip_v(6,lattice_maxNslip,lattice_Nstructure)); lattice_Sslip_v = 0.0_pReal
- allocate(lattice_sd(3,lattice_maxNslip,lattice_Nstructure)); lattice_sd = 0.0_pReal
+  allocate(lattice_sd(3,lattice_maxNslip,lattice_Nstructure)); lattice_sd = 0.0_pReal
- allocate(lattice_st(3,lattice_maxNslip,lattice_Nstructure)); lattice_st = 0.0_pReal
+  allocate(lattice_st(3,lattice_maxNslip,lattice_Nstructure)); lattice_st = 0.0_pReal
- allocate(lattice_sn(3,lattice_maxNslip,lattice_Nstructure)); lattice_sn = 0.0_pReal
+  allocate(lattice_sn(3,lattice_maxNslip,lattice_Nstructure)); lattice_sn = 0.0_pReal
 
- allocate(lattice_Qtwin(3,3,lattice_maxNtwin,lattice_Nstructure)); lattice_Qtwin = 0.0_pReal
+  allocate(lattice_Qtwin(3,3,lattice_maxNtwin,lattice_Nstructure)); lattice_Qtwin = 0.0_pReal
- allocate(lattice_Stwin(3,3,lattice_maxNtwin,lattice_Nstructure)); lattice_Stwin = 0.0_pReal
+  allocate(lattice_Stwin(3,3,lattice_maxNtwin,lattice_Nstructure)); lattice_Stwin = 0.0_pReal
- allocate(lattice_Stwin_v(6,lattice_maxNtwin,lattice_Nstructure)); lattice_Stwin_v = 0.0_pReal
+  allocate(lattice_Stwin_v(6,lattice_maxNtwin,lattice_Nstructure)); lattice_Stwin_v = 0.0_pReal
- allocate(lattice_td(3,lattice_maxNtwin,lattice_Nstructure)); lattice_td = 0.0_pReal
+  allocate(lattice_td(3,lattice_maxNtwin,lattice_Nstructure)); lattice_td = 0.0_pReal
- allocate(lattice_tt(3,lattice_maxNtwin,lattice_Nstructure)); lattice_tt = 0.0_pReal
+  allocate(lattice_tt(3,lattice_maxNtwin,lattice_Nstructure)); lattice_tt = 0.0_pReal
- allocate(lattice_tn(3,lattice_maxNtwin,lattice_Nstructure)); lattice_tn = 0.0_pReal
+  allocate(lattice_tn(3,lattice_maxNtwin,lattice_Nstructure)); lattice_tn = 0.0_pReal
- allocate(lattice_shearTwin(lattice_maxNtwin,lattice_Nstructure)); lattice_shearTwin = 0.0_pReal
+  allocate(lattice_shearTwin(lattice_maxNtwin,lattice_Nstructure)); lattice_shearTwin = 0.0_pReal
 
- allocate(lattice_interactionSlipSlip(lattice_maxNslip,lattice_maxNslip,lattice_Nstructure)); lattice_interactionSlipSlip = 0_pInt
+  allocate(lattice_interactionSlipSlip(lattice_maxNslip,lattice_maxNslip,lattice_Nstructure)); lattice_interactionSlipSlip = 0_pInt
- allocate(lattice_interactionSlipTwin(lattice_maxNslip,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionSlipTwin = 0_pInt
+  allocate(lattice_interactionSlipTwin(lattice_maxNslip,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionSlipTwin = 0_pInt
- allocate(lattice_interactionTwinTwin(lattice_maxNtwin,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionTwinTwin = 0_pInt
+  allocate(lattice_interactionTwinTwin(lattice_maxNtwin,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionTwinTwin = 0_pInt
- allocate(lattice_interactionTwinSlip(lattice_maxNtwin,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionTwinSlip = 0_pInt
+  allocate(lattice_interactionTwinSlip(lattice_maxNtwin,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionTwinSlip = 0_pInt
   
-end subroutine
+endsubroutine
 
 
 function lattice_initializeStructure(struct,CoverA)
@@ -758,7 +758,7 @@ function lattice_initializeStructure(struct,CoverA)
 
  lattice_initializeStructure = myStructure
 
-end function
+endfunction
 
 
 END MODULE

trunk/material.f90

@@ -102,7 +102,7 @@ subroutine material_init()
 
  call material_populateGrains()
 
-end subroutine
+endsubroutine
 
 
 !*********************************************************************
@@ -170,7 +170,7 @@ subroutine material_parseHomogenization(file,myPart)
 100 homogenization_maxNgrains = maxval(homogenization_Ngrains)
   return
 
- end subroutine
+ endsubroutine
 
 
 !*********************************************************************
@@ -240,7 +240,7 @@ subroutine material_parseMicrostructure(file,myPart)
 
 100 return
 
- end subroutine
+ endsubroutine
 
 
 !*********************************************************************
@@ -302,7 +302,7 @@ subroutine material_parsePhase(file,myPart)
 
 100 return
 
- end subroutine
+ endsubroutine
 
 
 !*********************************************************************
@@ -418,7 +418,7 @@ subroutine material_parseTexture(file,myPart)
 
 100 return
 
- end subroutine
+ endsubroutine
 
 
 !*********************************************************************
@@ -499,7 +499,7 @@ subroutine material_populateGrains()
            endif
          enddo
          NgrainsOfConstituent(t) = NgrainsOfConstituent(t) + sgn               ! change that by one
-       end do
+       enddo
 ! ----------------------------------------------------------------------------
        phaseOfGrain = 0_pInt
        orientationOfGrain = 0.0_pReal
@@ -563,7 +563,7 @@ subroutine material_populateGrains()
          endif
          
          grain = grain + NgrainsOfConstituent(i)                              ! advance microstructure grain index
-       end do  ! constituent
+       enddo  ! constituent
 
 ! ----------------------------------------------------------------------------
        do i=1,myNgrains-1                                                     ! walk thru grains
@@ -575,7 +575,7 @@ subroutine material_populateGrains()
          orientation             = orientationOfGrain(:,t)
          orientationOfGrain(:,t) = orientationOfGrain(:,i)
          orientationOfGrain(:,i) = orientation
-       end do
+       enddo
        !calc fraction after weighing with volumePerGrain
        !exchange in MC steps to improve result...
 
@@ -606,7 +606,7 @@ subroutine material_populateGrains()
  
  return
 
- end subroutine
+ endsubroutine
 
 
 END MODULE

trunk/mesh.f90

@@ -181,15 +181,16 @@
 !***********************************************************
 ! initialization 
 !***********************************************************
- SUBROUTINE mesh_init ()
+ subroutine mesh_init ()
 
  use prec, only: pInt
  use IO, only: IO_error,IO_open_InputFile
- use FEsolving, only: parallelExecution
+ use FEsolving, only: parallelExecution, FEsolving_execElem, FEsolving_execIP
  
  implicit none
  
  integer(pInt), parameter :: fileUnit = 222
+ integer(pInt) e
 
  mesh_Nelems          = 0_pInt
  mesh_NcpElems        = 0_pInt
@@ -228,14 +229,18 @@
    call IO_error(100) ! cannot open input file
  endif
  
- END SUBROUTINE
+ FEsolving_execElem = (/1,mesh_NcpElems/)
+ allocate(FEsolving_execIP(2,mesh_NcpElems)); FEsolving_execIP = 1_pInt
+ forall (e = 1:mesh_NcpElems) FEsolving_execIP(2,e) = FE_Nips(mesh_element(2,e))
+ 
+ endsubroutine
  
 
 
 !***********************************************************
 ! mapping of FE element types to internal representation
 !***********************************************************
- FUNCTION FE_mapElemtype(what)
+ function FE_mapElemtype(what)
 
  implicit none
  
@@ -259,9 +264,9 @@
        FE_mapElemtype = 7            ! Three-dimensional Arbitrarily Distorted qudratic hexahedral
 	case default 
 	   FE_mapElemtype = 0            ! unknown element --> should raise an error upstream..!
- end select
+ endselect
 
- END FUNCTION
+ endfunction
 
 
 
@@ -270,7 +275,7 @@
 !
 ! valid questions are 'elem', 'node'
 !***********************************************************
- FUNCTION mesh_FEasCP(what,id)
+ function mesh_FEasCP(what,id)
  
  use prec, only: pInt
  use IO, only: IO_lc
@@ -289,7 +294,7 @@
      lookupMap => mesh_mapFEtoCPnode
    case default
      return
- end select
+ endselect
  
  lower = 1_pInt
  upper = size(lookupMap,2)
@@ -313,17 +318,17 @@
    else
      mesh_FEasCP = lookupMap(2,center)
      exit
-   end if
+   endif
- end do
+ enddo
  return
  
- END FUNCTION
+ endfunction
 
 
 !***********************************************************
 ! find face-matching element of same type
 !!***********************************************************
- FUNCTION mesh_faceMatch(face,elem)
+ function mesh_faceMatch(face,elem)
 
  use prec, only: pInt
  implicit none
@@ -344,7 +349,7 @@
      minN = NsharedElems       ! remember min # shared elems
      lonelyNode = faceNode     ! remember most lonely node
    endif
- end do
+ enddo
 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 ?
@@ -358,13 +363,13 @@ candidate: do i=1,minN  ! iterate over lonelyNode's shared elements
        mesh_faceMatch = 0_pInt         ! set to "no match" (so far)
        cycle candidate                 ! next candidate elem
      endif
-   end do
+   enddo
    exit        ! surviving candidate
- end do candidate
+ enddo candidate
   
  return
  
- END FUNCTION
+ endfunction
 
  
 !********************************************************************
@@ -373,7 +378,7 @@ candidate: do i=1,minN  ! iterate over lonelyNode's shared elements
 ! assign globals:
 ! FE_nodesAtIP, FE_ipNeighbor, FE_subNodeParent, FE_subNodeOnIPFace
 !********************************************************************
- SUBROUTINE mesh_get_FEdata ()
+ subroutine mesh_get_FEdata ()
  
  use prec, only: pInt
  implicit none
@@ -1007,7 +1012,7 @@ candidate: do i=1,minN  ! iterate over lonelyNode's shared elements
  
  return
  
- END SUBROUTINE
+ endsubroutine
 
 !********************************************************************
 ! get count of elements, nodes, and cp elements in mesh
@@ -1016,7 +1021,7 @@ candidate: do i=1,minN  ! iterate over lonelyNode's shared elements
 ! assign globals:
 ! _Nelems, _Nnodes, _NcpElems
 !********************************************************************
- SUBROUTINE mesh_get_meshDimensions (unit)
+ subroutine mesh_get_meshDimensions (unit)
 
  use prec, only: pInt
  use IO
@@ -1046,19 +1051,19 @@ candidate: do i=1,minN  ! iterate over lonelyNode's shared elements
           case('element')                                ! Count the number of encountered element sets
             mesh_NelemSets=mesh_NelemSets+1
             mesh_maxNelemInSet = max(mesh_maxNelemInSet,IO_countContinousIntValues(unit))
-        end select
+        endselect
      case('hypoelastic')
        do i=1,3+hypoelasticTableStyle  ! Skip 3 or 4 lines
          read (unit,610,END=620) line
-       end do
+       enddo
        mesh_NcpElems = mesh_NcpElems + IO_countContinousIntValues(unit)
-   end select
+   endselect
 
- end do
+ enddo
 
 620 return
  
- END SUBROUTINE
+ endsubroutine
 
  
 !!********************************************************************
@@ -1068,7 +1073,7 @@ candidate: do i=1,minN  ! iterate over lonelyNode's shared elements
 ! assign globals:
 ! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsharedElems
 !********************************************************************
-SUBROUTINE mesh_get_nodeElemDimensions (unit)
+subroutine mesh_get_nodeElemDimensions (unit)
  
  use prec, only: pInt
  use IO
@@ -1107,25 +1112,25 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
            n = mesh_FEasCP('node',IO_IntValue (line,pos,j+2))
            if (all(node_seen /= n)) node_count(n) = node_count(n)+1
            node_seen(j) = n
-         end do
+         enddo
          call IO_skipChunks(unit,FE_NoriginalNodes(t)-(pos(1)-2))        ! read on if FE_Nnodes exceeds node count present on current line
-	   end if
+	   endif
-     end do
+     enddo
      exit
-   end if
+   endif
- end do
+ enddo
  
 630 mesh_maxNsharedElems = maxval(node_count)
  
  return
- END SUBROUTINE
+ endsubroutine
 
 !********************************************************************
 ! Build element set mapping 
 !
 ! allocate globals: mesh_nameElemSet, mesh_mapElemSet
 !********************************************************************
- SUBROUTINE mesh_build_elemSetMapping (unit)
+ subroutine mesh_build_elemSetMapping (unit)
 
  use prec, only: pInt
  use IO
@@ -1151,10 +1156,10 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
       elem_set = elem_set+1
       mesh_nameElemSet(elem_set) = IO_stringValue(line,pos,4)
       mesh_mapElemSet(:,elem_set) = IO_continousIntValues(unit,mesh_maxNelemInSet,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets)
-   end if
+   endif
- end do
+ enddo
 640 return
- END SUBROUTINE
+ endsubroutine
 
  
 !********************************************************************
@@ -1163,7 +1168,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
 ! allocate globals:
 ! _mapFEtoCPnode
 !********************************************************************
- SUBROUTINE mesh_build_nodeMapping (unit)
+ subroutine mesh_build_nodeMapping (unit)
 
  use prec, only: pInt
  use math, only: qsort
@@ -1191,15 +1196,15 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
        read (unit,610,END=650) line
        mesh_mapFEtoCPnode(1,i) = IO_fixedIntValue (line,(/0,10/),1)
        mesh_mapFEtoCPnode(2,i) = i
-     end do
+     enddo
      exit
-   end if
+   endif
- end do
+ enddo
 
 650 call qsort(mesh_mapFEtoCPnode,1,size(mesh_mapFEtoCPnode,2))
 
  return
- END SUBROUTINE
+ endsubroutine
 
 
 !********************************************************************
@@ -1208,7 +1213,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
 ! allocate globals:
 ! _mapFEtoCPelem
 !********************************************************************
- SUBROUTINE mesh_build_elemMapping (unit)
+ subroutine mesh_build_elemMapping (unit)
 
  use prec, only: pInt
  use math, only: qsort
@@ -1234,20 +1239,20 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
    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=660) line 
-     end do
+     enddo
      contInts = IO_continousIntValues(unit,mesh_NcpElems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets)
      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
+   endif
- end do
+ enddo
 
 660 call qsort(mesh_mapFEtoCPelem,1,size(mesh_mapFEtoCPelem,2))  ! should be mesh_NcpElems
 
  return
- END SUBROUTINE
+ endsubroutine
 
 
 !********************************************************************
@@ -1256,7 +1261,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
 ! allocate globals:
 ! _node
 !********************************************************************
- SUBROUTINE mesh_build_nodes (unit)
+ subroutine mesh_build_nodes (unit)
 
  use prec, only: pInt
  use IO
@@ -1283,15 +1288,15 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
        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
+       enddo
-     end do
+     enddo
      exit
-   end if
+   endif
- end do
+ enddo
 
 670 return
 
- END SUBROUTINE
+ endsubroutine
 
 
 !********************************************************************
@@ -1300,7 +1305,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
 ! allocate globals:
 ! _element
 !********************************************************************
- SUBROUTINE mesh_build_elements (unit)
+ subroutine mesh_build_elements (unit)
 
  use prec, only: pInt
  use IO
@@ -1334,8 +1339,8 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
          forall (j=1:FE_Nnodes(mesh_element(2,e))) &
            mesh_element(j+4,e) = IO_IntValue (line,pos,j+2) ! copy FE ids of nodes
          call IO_skipChunks(unit,FE_NoriginalNodes(mesh_element(2,e))-(pos(1)-2))        ! read on if FE_Nnodes exceeds node count present on current line
-	   end if
+	   endif
-     end do
+     enddo
 
      exit
    endif
@@ -1378,7 +1383,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
 
 620 return
 
- END SUBROUTINE
+ endsubroutine
 
  
 !********************************************************************
@@ -1387,7 +1392,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
 ! allocate globals:
 ! _sharedElem
 !********************************************************************
- SUBROUTINE mesh_build_sharedElems (unit)
+ subroutine mesh_build_sharedElems (unit)
 
  use prec, only: pInt
  use IO
@@ -1422,19 +1427,19 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
            if (all(node_seen /= n)) then
 		     mesh_sharedElem(1,n) = mesh_sharedElem(1,n) + 1
 		     mesh_sharedElem(1+mesh_sharedElem(1,n),n) = e
-           end if
+           endif
            node_seen(j) = n
 		 enddo
          call IO_skipChunks(unit,FE_NoriginalNodes(mesh_element(2,e))-(pos(1)-2))        ! read on if FE_Nnodes exceeds node count present on current line
-	   end if
+	   endif
-     end do
+     enddo
      exit
-   end if
+   endif
- end do
+ enddo
 
 620 return
 
- END SUBROUTINE
+ endsubroutine
  
 
 !***********************************************************
@@ -1443,7 +1448,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
 ! allocate globals
 ! _ipNeighborhood
 !***********************************************************
- SUBROUTINE mesh_build_ipNeighborhood()
+ subroutine mesh_build_ipNeighborhood()
  
  use prec, only: pInt
  implicit none
@@ -1508,7 +1513,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
  
  return
 
- END SUBROUTINE
+ endsubroutine
 
 
 
@@ -1518,7 +1523,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
 ! allocate globals
 ! _subNodeCoord
 !***********************************************************
- SUBROUTINE mesh_build_subNodeCoords()
+ subroutine mesh_build_subNodeCoords()
  
  use prec, only: pInt,pReal
  implicit none
@@ -1545,7 +1550,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
 
  return
  
- END SUBROUTINE
+ endsubroutine
 
 
 !***********************************************************
@@ -1554,7 +1559,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
 ! allocate globals
 ! _ipVolume
 !***********************************************************
- SUBROUTINE mesh_build_ipVolumes()
+ subroutine mesh_build_ipVolumes()
  
  use prec, only: pInt
  use math, only: math_volTetrahedron
@@ -1578,8 +1583,8 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
        do n = 1,FE_NipFaceNodes          ! loop over nodes on interface
          gravityNode(FE_subNodeOnIPFace(n,f,i,t)) = 1
          gravityNodePos(:,FE_subNodeOnIPFace(n,f,i,t)) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e)
-       end do
+       enddo
-     end do
+     enddo
      
      do j = 1,mesh_maxNnodes+mesh_maxNsubNodes-1        ! walk through entire flagList except last
        if (gravityNode(j) > 0_pInt) then                ! valid node index
@@ -1588,10 +1593,10 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
              gravityNode(j) = 0_pInt                    ! delete first instance
              gravityNodePos(:,j) = 0.0_pReal
              exit                                       ! continue with next suspect
-           end if
+           endif
-         end do
+         enddo
-       end if
+       endif
-     end do
+     enddo
      centerOfGravity = sum(gravityNodePos,2)/count(gravityNode > 0)
 
      do f = 1,FE_NipNeighbors(t)         ! loop over interfaces of IP and add tetrahedra which connect to CoG
@@ -1602,13 +1607,13 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
                                            nPos(:,1+mod(n+j-0,FE_NipFaceNodes)), &
                                            centerOfGravity)
        mesh_ipVolume(i,e) = mesh_ipVolume(i,e) + sum(volume)    ! add contribution from this interface
-     end do
+     enddo
      mesh_ipVolume(i,e) = mesh_ipVolume(i,e) / FE_NipFaceNodes  ! renormalize with interfaceNodeNum due to loop over them
-   end do
+   enddo
- end do
+ enddo
  return
  
- END SUBROUTINE
+ endsubroutine
 
 
 !***********************************************************
@@ -1617,7 +1622,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
 ! allocate globals
 ! _ipArea, _ipAreaNormal
 !***********************************************************
- SUBROUTINE mesh_build_ipAreas()
+ subroutine mesh_build_ipAreas()
  
  use prec, only: pInt,pReal
  use math
@@ -1652,7 +1657,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
  enddo
  return
  
- END SUBROUTINE
+ endsubroutine
 
 
 !***********************************************************
@@ -1660,7 +1665,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
 ! to the output file
 ! 
 !***********************************************************
- SUBROUTINE mesh_tell_statistics()
+ subroutine mesh_tell_statistics()
 
  use prec, only: pInt
  use math, only: math_range
@@ -1707,9 +1712,9 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
      write (6,"(i5,x,i5,x,e15.8)") e,i,mesh_IPvolume(i,e)
      do f = 1,FE_NipNeighbors(mesh_element(2,e))
        write (6,"(i33,x,e15.8,x,3(f6.3,x))") f,mesh_ipArea(f,i,e),mesh_ipAreaNormal(:,f,i,e)
-     end do
+     enddo
-   end do
+   enddo
- end do
+ enddo
  !write (6,*)
  !write (6,*) "Input Parser: SUBNODE COORDINATES"
  !write (6,*)
@@ -1723,10 +1728,10 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
  !                                                           mesh_subNodeCoord(1,FE_subNodeOnIPFace(n,f,i,t),e),&
  !                                                           mesh_subNodeCoord(2,FE_subNodeOnIPFace(n,f,i,t),e),&
  !                                                           mesh_subNodeCoord(3,FE_subNodeOnIPFace(n,f,i,t),e)
- !      end do
+ !      enddo
- !    end do
+ !    enddo
- !  end do
+ !  enddo
- !end do
+ !enddo
  write (6,*)
  write (6,*)
  write (6,*) "Input Parser: STATISTICS"
@@ -1757,7 +1762,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
 
  return
  
- END SUBROUTINE
+ endsubroutine
  
  
  END MODULE mesh

trunk/mpie_cpfem_marc.f90

@@ -16,6 +16,8 @@
 !             - set statevariable 3 to index of microstructure
 !             - make sure the file "material.config" exists in the working
 !               directory
+!             - make sure the file "numerics.config" exists in the working 
+!               directory
 !             - use nonsymmetric option for solver (e.g. direct 
 !               profile or multifrontal sparse, the latter seems
 !               to be faster!)
@@ -33,9 +35,10 @@
 !********************************************************************
 !
  include "prec.f90"             ! uses nothing else
- include "debug.f90"            ! uses prec
  include "math.f90"             ! uses prec
- include "IO.f90"               ! uses prec, debug, math
+ include "IO.f90"               ! uses prec, math
+ include "numerics.f90"         ! uses prec, IO
+ 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
@@ -44,10 +47,10 @@
  include "constitutive_j2.f90"                   ! uses prec, math, IO, lattice, material, debug
  include "constitutive_dislobased.f90"           ! uses prec, math, IO, lattice, material, debug
  include "constitutive.f90"     ! uses prec, IO, math, lattice, mesh, debug
- include "crystallite.f90"      ! uses
+ include "crystallite.f90"      ! uses prec, math, IO, numerics 
- include "homogenization_isostrain.f90"          ! uses
+ include "homogenization_isostrain.f90"          ! uses prec, math, IO, 
- include "homogenization.f90"   ! uses
+ include "homogenization.f90"   ! uses prec, math, IO, numerics
- include "CPFEM.f90"            ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite
+ include "CPFEM.f90"            ! uses prec, math, IO, numerics, debug, FEsolving, mesh, lattice, constitutive, crystallite
 
 
 !********************************************************************
@@ -123,11 +126,21 @@ subroutine hypela2(&
      ifu &        ! set to 1 if stretch has been calculated
    )
 
- use prec, only: pReal,pInt, iJacoStiffness
+ use prec, only:      pReal, &
- use FEsolving
+                      pInt
- use CPFEM, only: CPFEM_general
+ use FEsolving, only: cycleCounter, &
- use math, only: invnrmMandel
+                      theInc, &
- use debug
+                      theCycle, &
+                      theLovl, &
+                      theTime, &
+                      lastIncConverged, &
+                      outdatedByNewInc, &
+                      outdatedFFN1, &
+                      symmetricSolver
+ use CPFEM, only:     CPFEM_general
+ use math, only:      invnrmMandel
+ use debug, only:     debug_info, &
+                      debug_reset
  implicit none
  
 !     ** Start of generated type statements **
@@ -147,7 +160,7 @@ subroutine hypela2(&
  include "concom%%MARCVERSION%%"     ! concom is needed for inc, subinc, ncycle, lovl
  include "creeps%%MARCVERSION%%"     ! creeps is needed for timinc (time increment)
 
- integer(pInt) computationMode,i
+ integer(pInt)        computationMode, i
 
  if (inc == 0) then
    cycleCounter = 4
@@ -186,7 +199,7 @@ subroutine hypela2(&
  theCycle = ncycle                                  ! record current cycle count
  theLovl  = lovl                                    ! record current lovl
 
- call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,s,mod(cycleCounter-4,4_pInt*iJacoStiffness)==0,d,ngens)
+ call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,s,d,ngens)
 
 !     Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13
 !     Marc:   11, 22, 33, 12, 23, 13

trunk/numerics.config

@@ -0,0 +1,18 @@
+
+relevantStrain			1.0e-7
+iJacoStiffness			1
+iJacoLpresiduum			1
+pert_Fg				1.0e-6
+nHomog				10
+nCryst				20
+nState				10
+nStress				40
+subStepMin			1.0e-3
+rTol_crystalliteState		1.0e-6
+rTol_crystalliteStress		1.0e-6
+aTol_crystalliteStress		1.0e-8
+
+resToler			1.0e-4
+resAbsol			1.0e+2
+resBound			1.0e+1
+NRiterMax			24

trunk/numerics.f90

@@ -0,0 +1,169 @@
+!##############################################################
+MODULE numerics
+!##############################################################
+
+use prec, only: pInt, pReal
+implicit none
+
+character(len=64), parameter :: numerics_configFile = 'numerics.config' ! name of configuration file
+integer(pInt)                   iJacoStiffness, &                       ! frequency of stiffness update
+                                iJacoLpresiduum, &                      ! frequency of Jacobian update of residuum in Lp
+                                nHomog, &                               ! homogenization loop limit
+                                nCryst, &                               ! crystallite loop limit (only for debugging info, real loop limit is "subStepMin")
+                                nState, &                               ! state loop limit
+                                nStress, &                              ! stress loop limit
+                                NRiterMax                               ! maximum number of GIA iteration
+real(pReal)                     relevantStrain, &                       ! strain increment considered significant
+                                pert_Fg, &                              ! strain perturbation for FEM Jacobi
+                                subStepMin, &                           ! minimum (relative) size of sub-step allowed during cutback in crystallite
+                                rTol_crystalliteState, &                ! relative tolerance in crystallite state loop 
+                                rTol_crystalliteStress, &               ! relative tolerance in crystallite stress loop
+                                aTol_crystalliteStress, &               ! absolute tolerance in crystallite stress loop
+                                resToler, &                             ! relative tolerance of residual in GIA iteration
+                                resAbsol, &                             ! absolute tolerance of residual in GIA iteration (corresponds to ~1 Pa)
+                                resBound                                ! relative maximum value (upper bound) for GIA residual
+                                
+
+CONTAINS
+ 
+!*******************************************
+!    initialization subroutine
+!*******************************************
+subroutine numerics_init()
+  
+  !*** variables and functions from other modules ***!
+  use prec, only:                             pInt, & 
+                                              pReal  
+  use IO, only:                               IO_error, &
+                                              IO_open_file, &
+                                              IO_isBlank, &
+                                              IO_stringPos, &
+                                              IO_stringValue, &
+                                              IO_lc, &
+                                              IO_floatValue, &
+                                              IO_intValue
+  
+  implicit none
+
+  !*** input variables ***!
+  
+  !*** output variables ***!
+  
+  !*** local variables ***!
+  integer(pInt), parameter ::                 fileunit = 300  
+  integer(pInt), parameter ::                 maxNchunks = 2
+  integer(pInt), dimension(1+2*maxNchunks) :: positions
+  character(len=64)                           tag
+  character(len=1024)                         line
+  
+  !*** global variables ***!
+  ! relevantStrain
+  ! iJacoStiffness
+  ! iJacoLpresiduum
+  ! pert_Fg
+  ! nHomog
+  ! nCryst
+  ! nState
+  ! nStress
+  ! subStepMin
+  ! rTol_crystalliteState
+  ! rTol_crystalliteStress
+  ! aTol_crystalliteStress
+  ! resToler
+  ! resAbsol
+  ! resBound
+  ! NRiterMax
+  
+  ! initialize all values to zero
+  relevantStrain          = 0.0_pReal
+  iJacoStiffness          = 0_pInt
+  iJacoLpresiduum         = 0_pInt
+  pert_Fg                 = 0.0_pReal
+  nHomog                  = 0_pInt
+  nCryst                  = 0_pInt
+  nState                  = 0_pInt
+  nStress                 = 0_pInt
+  subStepMin              = 0.0_pReal
+  rTol_crystalliteState   = 0.0_pReal
+  rTol_crystalliteStress  = 0.0_pReal
+  aTol_crystalliteStress  = 0.0_pReal
+  resToler                = 0.0_pReal
+  resAbsol                = 0.0_pReal
+  resBound                = 0.0_pReal
+  NRiterMax               = 0_pInt
+
+  ! try to open the config file and call error if corrupt
+  if(.not. IO_open_file(fileunit,numerics_configFile)) call IO_error (100)
+ 
+  line = ''
+  ! read variables from config file
+  do
+    read(fileunit,'(a1024)',END=100) line
+    if (IO_isBlank(line)) cycle                           ! skip empty lines
+    positions = IO_stringPos(line,maxNchunks)
+    tag = IO_lc(IO_stringValue(line,positions,1))         ! extract key
+    select case(tag)
+      case ('relevantstrain')
+            relevantStrain = IO_floatValue(line,positions,2)
+      case ('ijacostiffness')
+            iJacoStiffness = IO_intValue(line,positions,2)
+      case ('ijacolpresiduum')
+            iJacoLpresiduum = IO_intValue(line,positions,2)
+      case ('pert_fg')
+            pert_Fg = IO_floatValue(line,positions,2)
+      case ('nhomog')
+            nHomog = IO_intValue(line,positions,2)
+      case ('ncryst')
+            nCryst = IO_intValue(line,positions,2)
+      case ('nstate')
+            nState = IO_intValue(line,positions,2)
+      case ('nstress')
+            nStress = IO_intValue(line,positions,2)
+      case ('substepmin')
+            subStepMin = IO_floatValue(line,positions,2)
+      case ('rtol_crystallitestate')
+            rTol_crystalliteState = IO_floatValue(line,positions,2)
+      case ('rtol_crystallitestress')
+            rTol_crystalliteStress = IO_floatValue(line,positions,2)
+      case ('atol_crystallitestress')
+            aTol_crystalliteStress = IO_floatValue(line,positions,2)
+      case ('restoler')
+            resToler = IO_floatValue(line,positions,2)
+      case ('resabsol')
+            resAbsol = IO_floatValue(line,positions,2)
+      case ('resbound')
+            resBound = IO_floatValue(line,positions,2)
+      case ('nritermax')
+            NRiterMax = IO_intValue(line,positions,2)
+    end select
+  enddo
+
+  100 write(6,*)
+  ! sanity check  
+  if (relevantStrain <= 0.0_pReal)          call IO_error(260)
+  if (iJacoStiffness < 1_pInt)              call IO_error(261)
+  if (iJacoLpresiduum < 1_pInt)             call IO_error(262)
+  if (pert_Fg <= 0.0_pReal)                 call IO_error(263)
+  if (nHomog < 1_pInt)                      call IO_error(264)
+  if (nCryst < 1_pInt)                      call IO_error(265)
+  if (nState < 1_pInt)                      call IO_error(266)
+  if (nStress < 1_pInt)                     call IO_error(267)
+  if (subStepMin <= 0.0_pReal)              call IO_error(268)
+  if (rTol_crystalliteState <= 0.0_pReal)   call IO_error(269)
+  if (rTol_crystalliteStress <= 0.0_pReal)  call IO_error(270)
+  if (aTol_crystalliteStress <= 0.0_pReal)  call IO_error(271)
+  if (resToler <= 0.0_pReal)                call IO_error(272)
+  if (resAbsol <= 0.0_pReal)                call IO_error(273)
+  if (resBound <= 0.0_pReal)                call IO_error(274)
+  if (NRiterMax < 1_pInt)                   call IO_error(275)
+  
+  close(fileunit)
+  
+  write(6,*)
+  write(6,*) '<<<+-  numerics init  -+>>>'
+  write(6,*) '...done'
+  write(6,*)
+ 
+endsubroutine
+
+END MODULE numerics

trunk/prec.f90

@@ -10,30 +10,8 @@
  integer, parameter :: pInt  = selected_int_kind(9)           ! up to +- 1e9
  integer, parameter :: pLongInt  = 8                          ! should be 64bit
 
-
  type :: p_vec
      real(pReal), dimension(:), pointer :: p
  end type p_vec
 
-!    *** Strain increment considered significant ***
- real(pReal), parameter :: relevantStrain = 1.0e-7_pReal
-
-!    *** Numerical parameters ***
- integer(pInt), parameter :: iJacoStiffness   = 1_pInt        ! frequency of stiffness update
- integer(pInt), parameter :: iJacoLpresiduum  = 6_pInt        ! frequency of Jacobian update of residuum in Lp
- real(pReal),   parameter :: pert_Fg          = 1.0e-6_pReal  ! strain perturbation for FEM Jacobi
- integer(pInt), parameter :: nHomog           = 10_pInt       ! homogenization loop limit
- integer(pInt), parameter :: nCryst           = 20_pInt       ! crystallite loop limit (only for debugging info, real loop limit is "subStepMin")
- integer(pInt), parameter :: nState           = 10_pInt       ! state loop limit
- integer(pInt), parameter :: nStress          = 40_pInt       ! stress loop limit
- real(pReal),   parameter :: rTol_crystalliteState  = 1.0e-5_pReal ! relative tolerance in crystallite state loop 
- real(pReal),   parameter :: rTol_crystalliteStress = 1.0e-6_pReal ! relative tolerance in crystallite stress loop
- real(pReal),   parameter :: aTol_crystalliteStress = 1.0e-8_pReal ! absolute tolerance in crystallite stress loop
- real(pReal),   parameter :: subStepMin = 1.0e-3_pReal        ! minimum (relative) size of sub-step allowed during cutback in crystallite
-!
- real(pReal),   parameter :: resToler = 1.0e-4_pReal          ! relative tolerance of residual in GIA iteration
- real(pReal),   parameter :: resAbsol = 1.0e+2_pReal          ! absolute tolerance of residual in GIA iteration (corresponds to ~1 Pa)
- real(pReal),   parameter :: resBound = 1.0e+1_pReal          ! relative maximum value (upper bound) for GIA residual
- integer(pInt), parameter :: NRiterMax = 24_pInt              ! maximum number of GIA iteration
-
  END MODULE prec