introduction of a possibility to have homogeneous element (all ips in one element are identical, sort of reduced integration)

and bugs-fixing in crystallite.f90, homogenization_RGC.f90, numerics.f90.

code/crystallite.f90

@@ -300,7 +300,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
   real(pReal), dimension(3,3,3,3) ::                    dPdF_pos,dPdF_neg
 
   ! ------ initialize to starting condition ------
-  centralDifference = .true.
   
 !$OMP CRITICAL (write2out)
 !  write (6,*)

code/homogenization_RGC.f90

@@ -303,7 +303,7 @@ function homogenization_RGC_updateState(&
  allocate(resid(3*nIntFaceTot));  resid = 0.0_pReal
  allocate(tract(nIntFaceTot,3));  tract = 0.0_pReal
  allocate(relax(3*nIntFaceTot));  relax = state%p(1:3*nIntFaceTot)
- allocate(drelax(3*nIntFaceTot)); &
+ allocate(drelax(3*nIntFaceTot))
    drelax = state%p(1:3*nIntFaceTot) - state0%p(1:3*nIntFaceTot)
 !* Debugging the obtained state
  if (RGCdebug) then
@@ -525,7 +525,7 @@ function homogenization_RGC_updateState(&
  
 !* Construct the Jacobian matrix of the numerical viscosity tangent
  allocate(rmatrix(3*nIntFaceTot,3*nIntFaceTot)); rmatrix = 0.0_pReal
- forall (i=1,3*nIntFaceTot) &
+ forall (i=1:3*nIntFaceTot) &
    rmatrix(i,i) = viscModus_RGC*ratePower_RGC/dt*(abs(drelax(i)/dt))**(ratePower_RGC - 1.0_pReal)
 !* Debugging the global Jacobian matrix of numerical viscosity tangent
  if (RGCdebugJacobi) then

code/material.config

@@ -27,6 +27,7 @@ Ngrains	2
 #####################
 
 [Aluminum_Poly]
+/elementhomogeneous/       # put this flag to set ips identical in one element (something like reduced integration)
 (constituent)	phase 3	texture	1	fraction	1.0
 
 [Aluminum_001]
@@ -36,6 +37,7 @@ Ngrains	2
 (constituent)	phase 1	texture	1	fraction	1.0
 
 [DP_Steel]
+/elementhomogeneous/
 (constituent)	phase	1	texture	1	fraction	0.82
 (constituent)	phase	2	texture	1	fraction	0.18 
 

code/material.f90

@@ -49,6 +49,7 @@ integer(pInt),     dimension(:),       allocatable :: homogenization_Ngrains, &
                                                       texture_Nfiber                   ! number of Fiber components per texture
 logical,           dimension(:),       allocatable :: homogenization_active, &         !
                                                       microstructure_active, &         ! 
+                                                      microstructure_elemhomo, &       ! flag to indicate homogeneous microstructure distribution over element's IPs
                                                       phase_localConstitution          ! flags phases with local constitutive law
 integer(pInt),     dimension(:,:),     allocatable :: microstructure_phase, &          ! phase IDs of each microstructure
                                                       microstructure_texture           ! texture IDs of each microstructure
@@ -102,9 +103,9 @@ subroutine material_init()
    write (6,'(x,a32,x,a16,x,i4)') homogenization_name(i),homogenization_type(i),homogenization_Ngrains(i)
  enddo
  write (6,*)
- write (6,'(a32,x,a12)') 'microstructure                  ','constituents'
+ write (6,'(a32,x,a12,x,a13)') 'microstructure                  ','constituents','homogeneous'
  do i = 1,material_Nmicrostructure
-   write (6,'(a32,x,i4)') microstructure_name(i),microstructure_Nconstituents(i)
+   write (6,'(a32,4x,i4,8x,l)') microstructure_name(i),microstructure_Nconstituents(i),microstructure_elemhomo(i)
    if (microstructure_Nconstituents(i) > 0_pInt) then
      do j = 1,microstructure_Nconstituents(i)
        write (6,'(a1,x,a32,x,a32,x,f6.4)') '>',phase_name(microstructure_phase(j,i)),&
@@ -207,14 +208,17 @@ subroutine material_parseMicrostructure(file,myPart)
  
  Nsections = IO_countSections(file,myPart)
  material_Nmicrostructure = Nsections
- allocate(microstructure_name(Nsections));    microstructure_name = ''
+ allocate(microstructure_name(Nsections));     microstructure_name = ''
  allocate(microstructure_Nconstituents(Nsections))
  allocate(microstructure_active(Nsections))
+ allocate(microstructure_elemhomo(Nsections))
 
  forall (i = 1:Nsections) microstructure_active(i) = any(mesh_element(4,:) == i)    ! current microstructure used in model?
 
  microstructure_Nconstituents = IO_countTagInPart(file,myPart,'(constituent)',Nsections)
  microstructure_maxNconstituents = maxval(microstructure_Nconstituents)
+ microstructure_elemhomo = IO_spotTagInPart(file,myPart,'/elementhomogeneous/',Nsections)
+
  allocate(microstructure_phase   (microstructure_maxNconstituents,Nsections)); microstructure_phase    = 0_pInt
  allocate(microstructure_texture (microstructure_maxNconstituents,Nsections)); microstructure_texture  = 0_pInt
  allocate(microstructure_fraction(microstructure_maxNconstituents,Nsections)); microstructure_fraction = 0.0_pReal
@@ -475,7 +479,12 @@ subroutine material_populateGrains()
      call IO_error(130,e,0,0)
    if (micro < 1 .or. micro > material_Nmicrostructure) &   ! out of bounds
      call IO_error(140,e,0,0)
-   Ngrains(homog,micro) = Ngrains(homog,micro) + homogenization_Ngrains(homog) * FE_Nips(mesh_element(2,e))
+   if (microstructure_elemhomo(micro)) then
+     dGrains = homogenization_Ngrains(homog)
+   else
+     dGrains = homogenization_Ngrains(homog) * FE_Nips(mesh_element(2,e))
+   endif
+   Ngrains(homog,micro) = Ngrains(homog,micro) + dGrains
  enddo
 
  allocate(volumeOfGrain(maxval(Ngrains)))           ! reserve memory for maximum case
@@ -499,10 +508,15 @@ subroutine material_populateGrains()
        grain = 0_pInt                               ! microstructure grain index
        do e = 1,mesh_NcpElems                       ! check each element
          if (mesh_element(3,e) == homog .and. mesh_element(4,e) == micro) then  ! my combination of homog and micro
-           forall (i = 1:FE_Nips(mesh_element(2,e))) &                          ! loop over IPs
+           if (microstructure_elemhomo(micro)) then                             ! homogeneous distribution of grains over each element's IPs
-             volumeOfGrain(grain+(i-1)*dGrains+1:grain+i*dGrains) = &
+             volumeOfGrain(grain+1:grain+dGrains) = sum(mesh_ipVolume(1:FE_Nips(mesh_element(2,e)),e))/dGrains
-               mesh_ipVolume(i,e)/dGrains                                       ! assign IPvolume/Ngrains to all grains of IP
+             grain = grain + dGrains                                            ! wind forward by NgrainsPerIP
-           grain = grain + FE_Nips(mesh_element(2,e)) * dGrains                 ! wind forward by Nips*NgrainsPerIP
+           else
+             forall (i = 1:FE_Nips(mesh_element(2,e))) &                        ! loop over IPs
+               volumeOfGrain(grain+(i-1)*dGrains+1:grain+i*dGrains) = &
+                 mesh_ipVolume(i,e)/dGrains                                     ! assign IPvolume/Ngrains to all grains of IP
+             grain = grain + FE_Nips(mesh_element(2,e)) * dGrains               ! wind forward by Nips*NgrainsPerIP
+           endif
          endif
        enddo
 ! ----------------------------------------------------------------------------  divide myNgrains as best over constituents
@@ -604,15 +618,24 @@ subroutine material_populateGrains()
        grain = 0_pInt                                                           ! microstructure grain index
        do e = 1,mesh_NcpElems                                                   ! check each element
          if (mesh_element(3,e) == homog .and. mesh_element(4,e) == micro) then  ! my combination of homog and micro
-           forall (i = 1:FE_Nips(mesh_element(2,e)), g = 1:dGrains)             ! loop over IPs and grains
+           if (microstructure_elemhomo(micro)) then                             ! homogeneous distribution of grains over each element's IPs
-             material_volume(g,i,e) = volumeOfGrain(grain+(i-1)*dGrains+g)
+             forall (i = 1:FE_Nips(mesh_element(2,e)), g = 1:dGrains)           ! loop over IPs and grains
-             material_phase(g,i,e) = phaseOfGrain(grain+(i-1)*dGrains+g)
+               material_volume(g,i,e) = volumeOfGrain(grain+g)
-             material_EulerAngles(:,g,i,e) = orientationOfGrain(:,grain+(i-1)*dGrains+g)
+               material_phase(g,i,e) = phaseOfGrain(grain+g)
-           end forall
+               material_EulerAngles(:,g,i,e) = orientationOfGrain(:,grain+g)
+             end forall
+             grain = grain + dGrains                                            ! wind forward by NgrainsPerIP
+           else
+             forall (i = 1:FE_Nips(mesh_element(2,e)), g = 1:dGrains)           ! loop over IPs and grains
+               material_volume(g,i,e) = volumeOfGrain(grain+(i-1)*dGrains+g)
+               material_phase(g,i,e) = phaseOfGrain(grain+(i-1)*dGrains+g)
+               material_EulerAngles(:,g,i,e) = orientationOfGrain(:,grain+(i-1)*dGrains+g)
+             end forall
+             grain = grain + FE_Nips(mesh_element(2,e)) * dGrains               ! wind forward by Nips*NgrainsPerIP
+           endif
            ! write (6,*) e
            ! write (6,*) material_phase(:,:,e)
            ! write (6,*) material_EulerAngles(:,:,:,e)
-           grain = grain + FE_Nips(mesh_element(2,e)) * dGrains                 ! wind forward by Nips*NgrainsPerIP
          endif
        enddo
 

code/numerics.f90

@@ -109,7 +109,7 @@ subroutine numerics_init()
   pPert_RGC               = 1.0e-7
   xSmoo_RGC               = 1.0e-5
   ratePower_RGC           = 1.0e+0  ! Newton viscosity (linear model)
-  viscModul_RGC           = 0.0e+0  ! No viscosity is applied
+  viscModus_RGC           = 0.0e+0  ! No viscosity is applied
   maxdRelax_RGC           = 1.0e+0
   
 !* Random seeding parameters: added <<<updated 27.08.2009>>>
@@ -186,7 +186,7 @@ subroutine numerics_init()
         case ('viscosityrate_rgc')
               ratePower_RGC = IO_floatValue(line,positions,2)
         case ('viscositymodulus_rgc')
-              viscModul_RGC = IO_floatValue(line,positions,2)
+              viscModus_RGC = IO_floatValue(line,positions,2)
         case ('maxrelaxation_rgc')
               maxdRelax_RGC = IO_floatValue(line,positions,2)
 
@@ -238,7 +238,7 @@ subroutine numerics_init()
   write(6,'(a24,x,e8.1)') 'perturbPenalty_RGC:   ',pPert_RGC
   write(6,'(a24,x,e8.1)') 'relevantMismatch_RGC: ',xSmoo_RGC
   write(6,'(a24,x,e8.1)') 'viscosityrate_RGC:    ',ratePower_RGC
-  write(6,'(a24,x,e8.1)') 'viscositymodulus_RGC: ',viscModul_RGC
+  write(6,'(a24,x,e8.1)') 'viscositymodulus_RGC: ',viscModus_RGC
   write(6,'(a24,x,e8.1)') 'maxrelaxation_RGC:    ',maxdRelax_RGC
   write(6,*)
 
@@ -277,7 +277,7 @@ subroutine numerics_init()
   if (pPert_RGC <= 0.0_pReal)               call IO_error(276)   !! oops !!
   if (xSmoo_RGC <= 0.0_pReal)               call IO_error(277)
   if (ratePower_RGC < 0.0_pReal)            call IO_error(278)
-  if (viscModul_RGC < 0.0_pReal)            call IO_error(278)
+  if (viscModus_RGC < 0.0_pReal)            call IO_error(278)
   if (maxdRelax_RGC <= 0.0_pReal)           call IO_error(288)
  
   if (fixedSeed <= 0_pInt)                  write(6,'(a)') 'Random is random!'