FEsolfing.f90: - added flag for symmetric solver

- added subroutine to detect symmetric solver

mesh.f90: - added subroutine call in mesh_init to detect symmetric solver during input file parsing

mpie_cpfem_marc2005r3.f90
mpie_cpfem_marc2007r1.f90
mpie_cpfem_marc2007r1_sequential.f90: - resorted include order of other source files
                                      - symmetrize d in case a symmetric solver is used

constitutive_pheno.f90: - included code to output shear rates and shear activity as post results

CPFEM_GIA8.f90
CPFEM_Taylor.f90
CPFEM_Taylor_sequential.f90: - symmetrize H_bar
                             - generalized reference to CPFEM_results in call of SingleCrystallite

trunk/CPFEM_GIA8.f90

@@ -136,7 +136,7 @@
 !
  integer(pInt) CPFEM_en, CPFEM_in, cp_en, CPFEM_ngens, i,j,k,l,m,n, e
  real(pReal), dimension (3,3)        :: ffn,ffn1,Kirchhoff_bar
- real(pReal), dimension (3,3,3,3)    :: H_bar
+ real(pReal), dimension (3,3,3,3)    :: H_bar, H_bar_sym
  real(pReal), dimension(CPFEM_ngens) :: CPFEM_stress
  real(pReal), dimension(CPFEM_ngens,CPFEM_ngens) :: CPFEM_jaco
  real(pReal) Temperature,CPFEM_dt,J_inverse
@@ -208,6 +208,8 @@
                            math_I3(j,l)*CPFEM_ffn1_bar(i,m,CPFEM_in,cp_en)*CPFEM_PK1_bar(k,m,CPFEM_in,cp_en)) + &
                           0.5_pReal*(math_I3(i,k)*Kirchhoff_bar(j,l) + math_I3(j,l)*Kirchhoff_bar(i,k) + &
                                      math_I3(i,l)*Kirchhoff_bar(j,k) + math_I3(j,k)*Kirchhoff_bar(i,l))
+       forall(i=1:3,j=1:3,k=1:3,l=1:3) &
+          H_bar_sym(i,j,k,l)= 0.25_pReal*(H_bar(i,j,k,l) +H_bar(j,i,k,l) +H_bar(i,j,l,k) +H_bar(j,i,l,k))
        CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel3333to66(J_inverse*H_bar)
 !
     case (3)    ! collect and return odd result
@@ -312,7 +314,8 @@
 ! -------------- grain loop -----------------
    do grain = 1,texture_Ngrains(mesh_element(4,cp_en))
      call SingleCrystallite(msg,PK1(:,:,grain),dPdF(:,:,:,:,grain),&
-                      CPFEM_results(5:4+constitutive_Nresults(grain,CPFEM_in,cp_en),grain,CPFEM_in,cp_en),&
+                      CPFEM_results(CPFEM_Nresults+1:CPFEM_Nresults+constitutive_Nresults(grain,CPFEM_in,cp_en),&
+                                    grain,CPFEM_in,cp_en),&
                       Fp1(:,:,grain),Fe1(:,:,grain),state1(:,grain),&   ! output up to here
                       dTime,cp_en,CPFEM_in,grain,.true.,&
                       CPFEM_Temperature(CPFEM_in,cp_en),F1(:,:,grain),F0(:,:,grain),Fp0(:,:,grain),state0(:,grain))
@@ -500,7 +503,8 @@
        call GIA_RelaxedDeformation(F1,F1_bar,rx)
        do grain = 1,8
          call SingleCrystallite(msg,PK1(:,:,grain),dPdF(:,:,:,:,grain),&
-                      CPFEM_results(5:4+constitutive_Nresults(grain,CPFEM_in,cp_en),grain,CPFEM_in,cp_en),&
+                      CPFEM_results(CPFEM_Nresults+1:CPFEM_Nresults+constitutive_Nresults(grain,CPFEM_in,cp_en),&
+                                    grain,CPFEM_in,cp_en),&
                       Fp1(:,:,grain),Fe1(:,:,grain),state1(:,grain),&   ! output up to here
                       dTime,cp_en,CPFEM_in,grain,.true.,&
                       CPFEM_Temperature(CPFEM_in,cp_en),F1(:,:,grain),F0(:,:,grain),Fp0(:,:,grain),state0(:,grain))

trunk/CPFEM_Taylor.f90

@@ -124,7 +124,7 @@
 !
  integer(pInt) CPFEM_en, CPFEM_in, cp_en, CPFEM_ngens, i,j,k,l,m,n, e
  real(pReal), dimension (3,3)        :: ffn,ffn1,Kirchhoff_bar
- real(pReal), dimension (3,3,3,3)    :: H_bar
+ real(pReal), dimension (3,3,3,3)    :: H_bar, H_bar_sym
  real(pReal), dimension(CPFEM_ngens) :: CPFEM_stress
  real(pReal), dimension(CPFEM_ngens,CPFEM_ngens) :: CPFEM_jaco, odd_jaco
  real(pReal) Temperature,CPFEM_dt,J_inverse
@@ -204,7 +204,8 @@
                           math_I3(j,l)*CPFEM_ffn1_bar(i,m,CPFEM_in,cp_en)*CPFEM_PK1_bar(k,m,CPFEM_in,cp_en) + &
                           0.5_pReal*(math_I3(i,k)*Kirchhoff_bar(j,l) + math_I3(j,l)*Kirchhoff_bar(i,k) + &
                                      math_I3(i,l)*Kirchhoff_bar(j,k) + math_I3(j,k)*Kirchhoff_bar(i,l))
-! Do we have to symmetrize H_bar?
+         forall(i=1:3,j=1:3,k=1:3,l=1:3) &
+            H_bar_sym(i,j,k,l)= 0.25_pReal*(H_bar(i,j,k,l) +H_bar(j,i,k,l) +H_bar(i,j,l,k) +H_bar(j,i,l,k))
          CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel3333to66(J_inverse*H_bar)
 ! if (CPFEM_in==8 .and. cp_en==80) then
 !   do e=1,80
@@ -288,7 +289,8 @@
  do grain = 1,texture_Ngrains(mesh_element(4,cp_en))
    dPdF = dPdF_bar_old                                                 ! preguess consistent tangent of grain with avg
    call SingleCrystallite(msg,PK1,dPdF,&
-                      CPFEM_results(5:4+constitutive_Nresults(grain,CPFEM_in,cp_en),grain,CPFEM_in,cp_en),&
+                      CPFEM_results(CPFEM_Nresults+1:CPFEM_Nresults+constitutive_Nresults(grain,CPFEM_in,cp_en),&
+                                    grain,CPFEM_in,cp_en),&
                       CPFEM_Lp(:,:,grain,CPFEM_in,cp_en),&
                       CPFEM_Fp_new(:,:,grain,CPFEM_in,cp_en),Fe1,constitutive_state_new(:,grain,CPFEM_in,cp_en),&   ! output up to here
                       CPFEM_dt,cp_en,CPFEM_in,grain,updateJaco,&

trunk/CPFEM_Taylor_sequential.f90

@@ -125,7 +125,7 @@
 !
  integer(pInt) CPFEM_en, CPFEM_in, cp_en, CPFEM_ngens, i,j,k,l,m,n, e
  real(pReal), dimension (3,3)        :: ffn,ffn1,Kirchhoff_bar
- real(pReal), dimension (3,3,3,3)    :: H_bar
+ real(pReal), dimension (3,3,3,3)    :: H_bar, H_bar_sym
  real(pReal), dimension(CPFEM_ngens) :: CPFEM_stress
  real(pReal), dimension(CPFEM_ngens,CPFEM_ngens) :: CPFEM_jaco, odd_jaco
  real(pReal) Temperature,CPFEM_dt,J_inverse
@@ -186,7 +186,9 @@
                         math_I3(j,l)*CPFEM_ffn1_bar(i,m,CPFEM_in,cp_en)*CPFEM_PK1_bar(k,m,CPFEM_in,cp_en) + &
                         0.5_pReal*(math_I3(i,k)*Kirchhoff_bar(j,l) + math_I3(j,l)*Kirchhoff_bar(i,k) + &
                                    math_I3(i,l)*Kirchhoff_bar(j,k) + math_I3(j,k)*Kirchhoff_bar(i,l))
-     CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel3333to66(J_inverse*H_bar)
+     forall(i=1:3,j=1:3,k=1:3,l=1:3) &
+        H_bar_sym(i,j,k,l)= 0.25_pReal*(H_bar(i,j,k,l) +H_bar(j,i,k,l) +H_bar(i,j,l,k) +H_bar(j,i,l,k))
+     CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel3333to66(J_inverse*H_bar)
 
    case (3)    ! *** collect and return odd result ***
      CPFEM_Temperature(CPFEM_in,cp_en)  = Temperature
@@ -264,7 +266,8 @@
    endif
 
    call SingleCrystallite(msg,PK1,dPdF,&
-                      CPFEM_results(5:4+constitutive_Nresults(grain,CPFEM_in,cp_en),grain,CPFEM_in,cp_en),&
+                      CPFEM_results(CPFEM_Nresults+1:CPFEM_Nresults+constitutive_Nresults(grain,CPFEM_in,cp_en),&
+                                    grain,CPFEM_in,cp_en),&
                       CPFEM_Lp(:,:,grain,CPFEM_in,cp_en),&
                       CPFEM_Fp_new(:,:,grain,CPFEM_in,cp_en),Fe1,constitutive_state_new(:,grain,CPFEM_in,cp_en),&   ! output up to here
                       CPFEM_dt,cp_en,CPFEM_in,grain,updateJaco,&

trunk/FEsolving.f90

@@ -10,5 +10,42 @@
  integer(pInt) theInc,theCycle,theLovl
  real(pReal)   theTime
  logical :: lastIncConverged = .false.,outdatedByNewInc = .false., outdatedFFN1 = .false.
+ logical :: symmetricSolver = .false. 
+
+ CONTAINS
+
+!***********************************************************
+! determine wether a symmetric solver is used 
+!***********************************************************
+ subroutine FE_get_solverSymmetry(unit)
+ 
+ use prec, only: pInt
+ use IO
+ implicit none
+ 
+ integer(pInt) unit
+ integer(pInt), dimension (133) :: pos
+ character*300 line
+ 
+610 FORMAT(A300)
+ 
+ rewind(unit)
+ do
+   read (unit,610,END=630) line
+   pos = IO_stringPos(line,1)
+   if( IO_lc(IO_stringValue(line,pos,1)) == 'solver' ) then
+     read (unit,610,END=630) line  ! Garbage line
+     pos = IO_stringPos(line,2)  ! limit to 64 nodes max (plus ID, type)
+     if(IO_intValue(line,pos,2) /= 1_pInt) symmetricSolver = .true.
+!$OMP CRITICAL (write2out)
+     write (6,*)
+     write (6,*) 'Symmetric solver detected. d-Matrix will be symmetrized!'
+!$OMP END CRITICAL (write2out)
+   endif
+ enddo
+
+630 return
+
+ end subroutine
 
  END MODULE FEsolving

trunk/constitutive_pheno.f90

@@ -701,7 +701,7 @@ do e=1,mesh_NcpElems
 		    constitutive_MatVolFrac(g,i,e) = 1.0_pReal ! singular material (so far)
 		    constitutive_TexVolFrac(g,i,e) = texVolfrac(s)/multiplicity(texID)/Nsym(texID)
 		    constitutive_Nstatevars(g,i,e) = material_Nslip(matID) ! number of state variables (i.e. tau_c of each slip system)
-		    constitutive_Nresults(g,i,e)   = 0         ! number of constitutive results
+!		    constitutive_Nresults(g,i,e)   = 2*material_Nslip(matID) ! number of constitutive results (shears in this case)
 		    constitutive_EulerAngles(:,g,i,e) = Euler(:,s)    ! store initial orientation
             forall (l=1:constitutive_Nstatevars(g,i,e))  ! initialize state variables
                constitutive_state_old(l,g,i,e) = material_s0_slip(matID)
@@ -900,7 +900,7 @@ implicit none
 !* Definition of variables
 integer(pInt) ipc,ip,el
 integer(pInt) matID,i
-real(pReal) dt,Temperature,tau_slip
+real(pReal) dt,Temperature,tau_slip, active_rate
 real(pReal), dimension(6) :: Tstar_v
 real(pReal), dimension(constitutive_Nstatevars(ipc,ip,el)) :: state
 real(pReal), dimension(constitutive_Nresults(ipc,ip,el))   :: constitutive_post_results
@@ -910,11 +910,18 @@ matID = constitutive_matID(ipc,ip,el)
 
 if(constitutive_Nresults(ipc,ip,el)==0) return
 
+constitutive_post_results=0
 do i=1,material_Nslip(matID)
-   constitutive_post_results(i) = state(i)
+!do i=1,constitutive_Nresults(ipc,ip,el)
+!   constitutive_post_results(i) = state(i)
    tau_slip=dot_product(Tstar_v,lattice_Sslip_v(:,i,material_CrystalStructure(matID)))
-   constitutive_post_results(i+material_Nslip(matID)) = &
+!   constitutive_post_results(i+material_Nslip(matID)) = &
-     dt*material_gdot0_slip(matID)*(abs(tau_slip)/state(i))**material_n_slip(matID)*sign(1.0_pReal,tau_slip)
+   constitutive_post_results(i) = &
+     material_gdot0_slip(matID)*(abs(tau_slip)/state(i))**material_n_slip(matID)*sign(1.0_pReal,tau_slip)
+enddo
+active_rate = 0.1_pReal*MAXVAL(abs(constitutive_post_results))
+do i=1,material_Nslip(matID)
+    if(abs(constitutive_post_results(i)) > active_rate) constitutive_post_results(i+material_Nslip(matID))=1.0_pReal
 enddo
 return
 

trunk/mesh.f90

@@ -220,7 +220,8 @@
  SUBROUTINE mesh_init ()
 
  use prec, only: pInt
- use IO, only: IO_error,IO_open_InputFile
+ use IO, only: IO_error,IO_open_InputFile
+ use FEsolving, only: FE_get_solverSymmetry
  implicit none
  
  integer(pInt), parameter :: fileUnit = 222
@@ -237,7 +238,8 @@
 
 
 ! call to various subroutines to parse the stuff from the input file...
- if (IO_open_inputFile(fileUnit)) then
+ if (IO_open_inputFile(fileUnit)) then
+   call FE_get_solverSymmetry(fileUnit)
    call mesh_get_meshDimensions(fileUnit)
    call mesh_build_nodeMapping(fileUnit)
    call mesh_build_elemMapping(fileUnit)

trunk/mpie_cpfem_marc2005r3.f90

@@ -27,16 +27,16 @@
 !             - creeps: timinc
 !********************************************************************
 !
- include "prec.f90"
+ include "prec.f90"             ! uses nothing else
- include "FEsolving.f90"
+ include "debug.f90"            ! uses prec
- include "debug.f90"
+ include "math.f90"             ! uses prec
- include "math.f90"
+ include "IO.f90"               ! uses prec, debug, math
- include "IO.f90"
+ include "FEsolving.f90"        ! uses prec, IO
- include "mesh.f90"
+ include "mesh.f90"             ! uses prec, IO, math, FEsolving
- include "lattice.f90"
+ include "lattice.f90"          ! uses prec, math
- include "constitutive.f90"
+ include "constitutive.f90"     ! uses prec, IO, math, lattice, mesh, debug
- include "crystallite.f90"
+ include "crystallite.f90"      ! uses prec, debug, constitutive, mesh, math, IO
- include "CPFEM.f90"
+ include "CPFEM.f90"            ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite
 !
 
  SUBROUTINE hypela2(d,g,e,de,s,t,dt,ngens,n,nn,kcus,matus,ndi,&
@@ -236,7 +236,7 @@
 !     Marc:   11, 22, 33, 12, 23, 13
  forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*d(1:ngens,i)*invnrmMandel(1:ngens)
  s(1:ngens) = s(1:ngens)*invnrmMandel(1:ngens)
-
+ if(symmetricSolver) d(1:ngens,1:ngens) = 0.5_pReal*(d(1:ngens,1:ngens)+transpose(d(1:ngens,1:ngens)))
  return
  
  END SUBROUTINE

trunk/mpie_cpfem_marc2007r1.f90

@@ -27,16 +27,16 @@
 !             - creeps: timinc
 !********************************************************************
 !
- include "prec.f90"
+ include "prec.f90"             ! uses nothing else
- include "FEsolving.f90"
+ include "debug.f90"            ! uses prec
- include "debug.f90"
+ include "math.f90"             ! uses prec
- include "math.f90"
+ include "IO.f90"               ! uses prec, debug, math
- include "IO.f90"
+ include "FEsolving.f90"        ! uses prec, IO
- include "mesh.f90"
+ include "mesh.f90"             ! uses prec, IO, math, FEsolving
- include "lattice.f90"
+ include "lattice.f90"          ! uses prec, math
- include "constitutive.f90"
+ include "constitutive.f90"     ! uses prec, IO, math, lattice, mesh, debug
- include "crystallite.f90"
+ include "crystallite.f90"      ! uses prec, debug, constitutive, mesh, math, IO
- include "CPFEM.f90"
+ include "CPFEM.f90"            ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite
 !
 
  SUBROUTINE hypela2(d,g,e,de,s,t,dt,ngens,n,nn,kcus,matus,ndi,&
@@ -213,7 +213,7 @@
 !     Marc:   11, 22, 33, 12, 23, 13
  forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*d(1:ngens,i)*invnrmMandel(1:ngens)
  s(1:ngens) = s(1:ngens)*invnrmMandel(1:ngens)
-
+ if(symmetricSolver) d(1:ngens,1:ngens) = 0.5_pReal*(d(1:ngens,1:ngens)+transpose(d(1:ngens,1:ngens)))
  return
  
  END SUBROUTINE
@@ -253,7 +253,7 @@
 ! assign result variable
  v=CPFEM_results(mod(jpltcd-1_pInt, CPFEM_Nresults+constitutive_maxNresults)+1_pInt,&
                  (jpltcd-1_pInt)/(CPFEM_Nresults+constitutive_maxNresults)+1_pInt,&
-                 nn, mesh_FEasCP('elem', m))
+                 nn, mesh_FEasCP('elem', m))
  return
  END SUBROUTINE
 !

trunk/mpie_cpfem_marc2007r1_sequential.f90

@@ -27,16 +27,16 @@
 !             - creeps: timinc
 !********************************************************************
 !
- include "prec.f90"
+ include "prec.f90"             ! uses nothing else
- include "FEsolving.f90"
+ include "debug.f90"            ! uses prec
- include "debug.f90"
+ include "math.f90"             ! uses prec
- include "math.f90"
+ include "IO.f90"               ! uses prec, debug, math
- include "IO.f90"
+ include "FEsolving.f90"        ! uses prec, IO
- include "mesh.f90"
+ include "mesh.f90"             ! uses prec, IO, math, FEsolving
- include "lattice.f90"
+ include "lattice.f90"          ! uses prec, math
- include "constitutive.f90"
+ include "constitutive.f90"     ! uses prec, IO, math, lattice, mesh, debug
- include "crystallite.f90"
+ include "crystallite.f90"      ! uses prec, debug, constitutive, mesh, math, IO
- include "CPFEM_sequential.f90"
+ include "CPFEM_sequential.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite
 !
 
  SUBROUTINE hypela2(d,g,e,de,s,t,dt,ngens,n,nn,kcus,matus,ndi,&
@@ -207,7 +207,7 @@
 !     Marc:   11, 22, 33, 12, 23, 13
  forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*d(1:ngens,i)*invnrmMandel(1:ngens)
  s(1:ngens) = s(1:ngens)*invnrmMandel(1:ngens)
-
+ if(symmetricSolver) d(1:ngens,1:ngens) = 0.5_pReal*(d(1:ngens,1:ngens)+transpose(d(1:ngens,1:ngens)))
  return
  
  END SUBROUTINE