DAMASK_EICMD/code/material.f90

! Copyright 2011 Max-Planck-Institut für Eisenforschung GmbH
!
! This file is part of DAMASK,
! the Düsseldorf Advanced MAterial Simulation Kit.
!
! DAMASK is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! DAMASK is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with DAMASK. If not, see <http://www.gnu.org/licenses/>.
!
!--------------------------------------------------------------------------------------------------
!* $Id$
!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!! Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Parses material.config
!--------------------------------------------------------------------------------------------------
module material

 use prec,  only: pReal, &
                  pInt

 implicit none
 private
 character(len=64), parameter, public  :: &
   material_configFile         = 'material.config', &
   material_localFileExt       = 'materialConfig'
   
 character(len=32), parameter, public  :: &
   material_partHomogenization = 'homogenization', &
   material_partCrystallite    = 'crystallite', &
   material_partPhase          = 'phase'
 
 character(len=64), dimension(:), allocatable, public :: &
   phase_elasticity, &                                                                              !> elasticity of each phase  
   phase_plasticity, &                                                                              !> plasticity of each phase  
   phase_name, &                                                                                    !> name of each phase
   homogenization_name, &                                                                           !> name of each homogenization
   homogenization_type, &                                                                           !> type of each homogenization
   crystallite_name                                                                                 !> name of each crystallite setting

 integer(pInt), public :: &
   homogenization_maxNgrains, &                                                                     !> max number of grains in any USED homogenization
   material_Nphase, &                                                                               !> number of phases
   material_Nhomogenization, &                                                                      !> number of homogenizations
   material_Nmicrostructure, &                                                                      !> number of microstructures
   material_Ncrystallite                                                                            !> number of crystallite settings
 
 integer(pInt), dimension(:), allocatable, public :: &
   homogenization_Ngrains, &                                                                        !> number of grains in each homogenization
   homogenization_Noutput, &                                                                        !> number of '(output)' items per homogenization
   phase_Noutput, &                                                                                 !> number of '(output)' items per phase
   phase_elasticityInstance, &                                                                      !> instance of particular elasticity of each phase
   phase_plasticityInstance, &                                                                      !> instance of particular plasticity of each phase
   crystallite_Noutput, &                                                                           !> number of '(output)' items per crystallite setting
   homogenization_typeInstance, &                                                                   !> instance of particular type of each homogenization
   microstructure_crystallite                                                                       !> crystallite setting ID of each microstructure

 integer(pInt), dimension(:,:,:), allocatable, public :: &
   material_phase, &                                                                                !> phase   (index) of each grain,IP,element
   material_texture                                                                                 !> texture (index) of each grain,IP,element
 
 real(pReal), dimension(:,:,:,:), allocatable, public :: &
   material_EulerAngles                                                                             !> initial orientation of each grain,IP,element
 
 logical, dimension(:), allocatable, public :: &
   microstructure_active, & 
   microstructure_elemhomo, &                                                                       !> flag to indicate homogeneous microstructure distribution over element's IPs
   phase_localPlasticity                                                                            !> flags phases with local constitutive law


 character(len=32), parameter, private :: &
   material_partMicrostructure = 'microstructure', &
   material_partTexture        = 'texture'
   
 character(len=64), dimension(:), allocatable, private :: &
   microstructure_name, &                                                                           !> name of each microstructure
   texture_name                                                                                     !> name of each texture
     
 character(len=256), dimension(:), allocatable, private :: &
   texture_ODFfile                                                                                  !> name of each ODF file         

 integer(pInt), private :: &
   material_Ntexture, &                                                                             !> number of textures
   microstructure_maxNconstituents, &                                                               !> max number of constituents in any phase
   texture_maxNgauss, &                                                                             !> max number of Gauss components in any texture
   texture_maxNfiber                                                                                !> max number of Fiber components in any texture

 integer(pInt), dimension(:), allocatable, private :: &
   microstructure_Nconstituents, &                                                                  !> number of constituents in each microstructure
   texture_symmetry, &                                                                              !> number of symmetric orientations per texture
   texture_Ngauss, &                                                                                !> number of Gauss components per texture
   texture_Nfiber                                                                                   !> number of Fiber components per texture
 
 integer(pInt), dimension(:,:), allocatable, private :: &
   microstructure_phase, &                                                                          !> phase IDs of each microstructure
   microstructure_texture                                                                           !> texture IDs of each microstructure
 
 real(pReal), dimension(:,:), allocatable, private :: &
   microstructure_fraction                                                                          !> vol fraction of each constituent in microstructure
 
 real(pReal), dimension(:,:,:), allocatable :: &
   material_volume, &                                                                               !> volume of each grain,IP,element
   texture_Gauss, &                                                                                 !> data of each Gauss component
   texture_Fiber                                                                                    !> data of each Fiber component
 
 logical, dimension(:), allocatable, private :: &
   homogenization_active


 public  :: material_init  
 
contains

!*********************************************************************
subroutine material_init
!*********************************************************************
!*      Module initialization         *
!**************************************

 use, intrinsic :: iso_fortran_env                                ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
 use IO,    only: IO_error, &
                  IO_open_file, &
                  IO_open_jobFile_stat
 use debug, only: debug_what, &
                  debug_material, &
                  debug_levelBasic, &
                  debug_levelExtensive
 
 implicit none

 integer(pInt), parameter :: fileunit = 200_pInt
 integer(pInt)            :: i,j, myDebug
 
 myDebug = debug_what(debug_material)
 
 !$OMP CRITICAL (write2out)
   write(6,*)
   write(6,*) '<<<+-  material init  -+>>>'
   write(6,*) '$Id$'
#include "compilation_info.f90"
 !$OMP END CRITICAL (write2out)
 
 if (.not. IO_open_jobFile_stat(fileunit,material_localFileExt)) then        ! no local material configuration present...
   call IO_open_file(fileunit,material_configFile)                           ! ...open material.config file
 endif
 call material_parseHomogenization(fileunit,material_partHomogenization)
 if (iand(myDebug,debug_levelBasic) /= 0_pInt) then
   !$OMP CRITICAL (write2out)
   write (6,*) 'Homogenization parsed'
   !$OMP END CRITICAL (write2out)
 endif
 call material_parseMicrostructure(fileunit,material_partMicrostructure)
 if (iand(myDebug,debug_levelBasic) /= 0_pInt) then
   !$OMP CRITICAL (write2out)
   write (6,*) 'Microstructure parsed'
   !$OMP END CRITICAL (write2out)
 endif
 call material_parseCrystallite(fileunit,material_partCrystallite)
 if (iand(myDebug,debug_levelBasic) /= 0_pInt) then
   !$OMP CRITICAL (write2out)
   write (6,*) 'Crystallite parsed'
   !$OMP END CRITICAL (write2out)
 endif
 call material_parseTexture(fileunit,material_partTexture)
 if (iand(myDebug,debug_levelBasic) /= 0_pInt) then
   !$OMP CRITICAL (write2out)
   write (6,*) 'Texture parsed'
   !$OMP END CRITICAL (write2out)
 endif
 call material_parsePhase(fileunit,material_partPhase)
 if (iand(myDebug,debug_levelBasic) /= 0_pInt) then
   !$OMP CRITICAL (write2out)
   write (6,*) 'Phase parsed'
   !$OMP END CRITICAL (write2out)
 endif
 close(fileunit)

 do i = 1_pInt,material_Nmicrostructure
   if (microstructure_crystallite(i) < 1_pInt .or. &
       microstructure_crystallite(i) > material_Ncrystallite) call IO_error(150_pInt,i)
   if (minval(microstructure_phase(1:microstructure_Nconstituents(i),i)) < 1_pInt .or. &
       maxval(microstructure_phase(1:microstructure_Nconstituents(i),i)) > material_Nphase) call IO_error(151_pInt,i)
   if (minval(microstructure_texture(1:microstructure_Nconstituents(i),i)) < 1_pInt .or. &
       maxval(microstructure_texture(1:microstructure_Nconstituents(i),i)) > material_Ntexture) call IO_error(152_pInt,i)
   if (abs(sum(microstructure_fraction(:,i)) - 1.0_pReal) >= 1.0e-10_pReal) then
     if (iand(myDebug,debug_levelExtensive) /= 0_pInt) then
       !$OMP CRITICAL (write2out)
         write(6,*)'sum of microstructure fraction = ',sum(microstructure_fraction(:,i))
       !$OMP END CRITICAL (write2out)
     endif
     call IO_error(153_pInt,i)
   endif
 enddo
 if (iand(myDebug,debug_levelExtensive) /= 0_pInt) then
   !$OMP CRITICAL (write2out)
     write (6,*)
     write (6,*) 'MATERIAL configuration'
     write (6,*)
     write (6,'(a32,1x,a16,1x,a6)') 'homogenization                  ','type            ','grains'
     do i = 1_pInt,material_Nhomogenization
       write (6,'(1x,a32,1x,a16,1x,i4)') homogenization_name(i),homogenization_type(i),homogenization_Ngrains(i)
     enddo
     write (6,*)
     write (6,'(a32,1x,a11,1x,a12,1x,a13)') 'microstructure                  ','crystallite','constituents','homogeneous'
     do i = 1_pInt,material_Nmicrostructure
       write (6,'(a32,4x,i4,8x,i4,8x,l1)') microstructure_name(i), &
                                           microstructure_crystallite(i), &
                                           microstructure_Nconstituents(i), &
                                           microstructure_elemhomo(i)
       if (microstructure_Nconstituents(i) > 0_pInt) then
         do j = 1_pInt,microstructure_Nconstituents(i)
           write (6,'(a1,1x,a32,1x,a32,1x,f7.4)') '>',phase_name(microstructure_phase(j,i)),&
                                                      texture_name(microstructure_texture(j,i)),&
                                                      microstructure_fraction(j,i)
         enddo
         write (6,*)
       endif
     enddo
   !$OMP END CRITICAL (write2out)
 endif
 
 call material_populateGrains

end subroutine material_init


!*********************************************************************
subroutine material_parseHomogenization(myFile,myPart)
!*********************************************************************

 use IO
 use mesh, only: mesh_element
 
 implicit none
 character(len=*), intent(in) :: myPart
 integer(pInt),    intent(in) :: myFile
 
 integer(pInt),     parameter :: maxNchunks = 2_pInt
 
 integer(pInt), dimension(1+2*maxNchunks) :: positions
 integer(pInt) Nsections, section, s
 character(len=64)   :: tag
 character(len=1024) ::line
 
 Nsections = IO_countSections(myFile,myPart)
 material_Nhomogenization = Nsections
 if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart)
 
 allocate(homogenization_name(Nsections));    homogenization_name = ''
 allocate(homogenization_type(Nsections));    homogenization_type = ''
 allocate(homogenization_typeInstance(Nsections));  homogenization_typeInstance = 0_pInt
 allocate(homogenization_Ngrains(Nsections)); homogenization_Ngrains = 0_pInt
 allocate(homogenization_Noutput(Nsections)); homogenization_Noutput = 0_pInt
 allocate(homogenization_active(Nsections));  homogenization_active = .false.

 forall (s = 1_pInt:Nsections) homogenization_active(s) = any(mesh_element(3,:) == s)    ! current homogenization used in model? Homogenization view, maximum operations depend on maximum number of homog schemes
 homogenization_Noutput = IO_countTagInPart(myFile,myPart,'(output)',Nsections)
 
 rewind(myFile)
 line = ''
 section = 0_pInt
 
 do while (IO_lc(IO_getTag(line,'<','>')) /= myPart)      ! wind forward to myPart
   read(myFile,'(a1024)',END=100) line
 enddo

 do
   read(myFile,'(a1024)',END=100) line
   if (IO_isBlank(line)) cycle                            ! skip empty lines
   if (IO_getTag(line,'<','>') /= '') exit                ! stop at next part
   if (IO_getTag(line,'[',']') /= '') then                ! next section
     section = section + 1_pInt
     homogenization_name(section) = IO_getTag(line,'[',']')
   endif
   if (section > 0_pInt) then
     positions = IO_stringPos(line,maxNchunks)
     tag = IO_lc(IO_stringValue(line,positions,1_pInt))        ! extract key
     select case(tag)
       case ('type')
         homogenization_type(section) = IO_lc(IO_stringValue(line,positions,2_pInt))  ! adding: IO_lc function <<<updated 31.07.2009>>>
         do s = 1_pInt,section
           if (homogenization_type(s) == homogenization_type(section)) &
             homogenization_typeInstance(section) = homogenization_typeInstance(section) + 1_pInt  ! count instances
         enddo
       case ('ngrains')
         homogenization_Ngrains(section) = IO_intValue(line,positions,2_pInt)
     end select
   endif
 enddo

100 homogenization_maxNgrains = maxval(homogenization_Ngrains,homogenization_active)

 end subroutine material_parseHomogenization


!*********************************************************************
subroutine material_parseMicrostructure(myFile,myPart)
!*********************************************************************

 use IO
 use mesh, only: mesh_element, mesh_NcpElems
 
 implicit none
 character(len=*), intent(in) :: myPart
 integer(pInt),    intent(in) :: myFile
 
 integer(pInt), parameter :: maxNchunks = 7_pInt
 
 integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions
 integer(pInt) :: Nsections, section, constituent, e, i
 character(len=64)   :: tag
 character(len=1024) :: line

 Nsections = IO_countSections(myFile,myPart)
 material_Nmicrostructure = Nsections
 if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart)

 allocate(microstructure_name(Nsections));            microstructure_name = ''
 allocate(microstructure_crystallite(Nsections));     microstructure_crystallite = 0_pInt
 allocate(microstructure_Nconstituents(Nsections))
 allocate(microstructure_active(Nsections))
 allocate(microstructure_elemhomo(Nsections))

 forall (e = 1_pInt:mesh_NcpElems) microstructure_active(mesh_element(4,e)) = .true.       ! current microstructure used in model? Elementwise view, maximum N operations for N elements
  
 microstructure_Nconstituents = IO_countTagInPart(myFile,myPart,'(constituent)',Nsections)
 microstructure_maxNconstituents = maxval(microstructure_Nconstituents)
 microstructure_elemhomo = IO_spotTagInPart(myFile,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
 
 rewind(myFile)
 line = ''
 section = 0_pInt
 
 do while (IO_lc(IO_getTag(line,'<','>')) /= myPart)      ! wind forward to myPart
   read(myFile,'(a1024)',END=100) line
 enddo

 do
   read(myFile,'(a1024)',END=100) line
   if (IO_isBlank(line)) cycle                            ! skip empty lines
   if (IO_getTag(line,'<','>') /= '') exit                ! stop at next part
   if (IO_getTag(line,'[',']') /= '') then                ! next section
     section = section + 1_pInt
     constituent = 0_pInt
     microstructure_name(section) = IO_getTag(line,'[',']')
   endif
   if (section > 0_pInt) then
     positions = IO_stringPos(line,maxNchunks)
     tag = IO_lc(IO_stringValue(line,positions,1_pInt))        ! extract key
     select case(tag)
       case ('crystallite')
         microstructure_crystallite(section) = IO_intValue(line,positions,2_pInt)
       case ('(constituent)')
         constituent = constituent + 1_pInt
         do i=2_pInt,6_pInt,2_pInt
           tag = IO_lc(IO_stringValue(line,positions,i))
           select case (tag)
             case('phase')
               microstructure_phase(constituent,section) =    IO_intValue(line,positions,i+1_pInt)
             case('texture')
               microstructure_texture(constituent,section) =  IO_intValue(line,positions,i+1_pInt)
             case('fraction')
               microstructure_fraction(constituent,section) = IO_floatValue(line,positions,i+1_pInt)
           end select
         enddo
     end select
   endif
 enddo

100 end subroutine material_parseMicrostructure


!*********************************************************************
subroutine material_parseCrystallite(myFile,myPart)
!*********************************************************************

 use IO, only: IO_countSections, &
               IO_error, &
               IO_countTagInPart, &
               IO_getTag, &
               IO_lc, &
               IO_isBlank

 implicit none
 character(len=*), intent(in) :: myPart
 integer(pInt),    intent(in) :: myFile
 
 integer(pInt)       :: Nsections, &
                        section
 character(len=1024) :: line
 
 Nsections = IO_countSections(myFile,myPart)
 material_Ncrystallite = Nsections
 if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart)

 allocate(crystallite_name(Nsections));       crystallite_name = ''
 allocate(crystallite_Noutput(Nsections));    crystallite_Noutput = 0_pInt

 crystallite_Noutput = IO_countTagInPart(myFile,myPart,'(output)',Nsections)
 
 rewind(myFile)
 line = ''
 section = 0_pInt
 
 do while (IO_lc(IO_getTag(line,'<','>')) /= myPart)      ! wind forward to myPart
   read(myFile,'(a1024)',END=100) line
 enddo

 do
   read(myFile,'(a1024)',END=100) line
   if (IO_isBlank(line)) cycle                            ! skip empty lines
   if (IO_getTag(line,'<','>') /= '') exit                ! stop at next part
   if (IO_getTag(line,'[',']') /= '') then                ! next section
     section = section + 1_pInt
     crystallite_name(section) = IO_getTag(line,'[',']')
   endif
 enddo

100 end subroutine material_parseCrystallite


!*********************************************************************
subroutine material_parsePhase(myFile,myPart)
!*********************************************************************

 use IO
 
 implicit none
 character(len=*), intent(in) :: myPart
 integer(pInt),    intent(in) :: myFile
 
 integer(pInt), parameter :: maxNchunks = 2_pInt
 
 integer(pInt), dimension(1+2*maxNchunks) :: positions
 integer(pInt) Nsections, section, s
 character(len=64)   :: tag
 character(len=1024) :: line
 
 Nsections = IO_countSections(myFile,myPart)
 material_Nphase = Nsections
 if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart)

 allocate(phase_name(Nsections));          phase_name = ''
 allocate(phase_elasticity(Nsections));  phase_elasticity = ''
 allocate(phase_elasticityInstance(Nsections));  phase_elasticityInstance = 0_pInt
 allocate(phase_plasticity(Nsections));  phase_plasticity = ''
 allocate(phase_plasticityInstance(Nsections));  phase_plasticityInstance = 0_pInt
 allocate(phase_Noutput(Nsections))
 allocate(phase_localPlasticity(Nsections))

 phase_Noutput = IO_countTagInPart(myFile,myPart,'(output)',Nsections)
 phase_localPlasticity = .not. IO_spotTagInPart(myFile,myPart,'/nonlocal/',Nsections)
 
 rewind(myFile)
 line = ''
 section = 0_pInt
 
 do while (IO_lc(IO_getTag(line,'<','>')) /= myPart)      ! wind forward to myPart
   read(myFile,'(a1024)',END=100) line
 enddo

 do
   read(myFile,'(a1024)',END=100) line
   if (IO_isBlank(line)) cycle                            ! skip empty lines
   if (IO_getTag(line,'<','>') /= '') exit                ! stop at next part
   if (IO_getTag(line,'[',']') /= '') then                ! next section
     section = section + 1_pInt
     phase_name(section) = IO_getTag(line,'[',']')
   endif
   if (section > 0_pInt) then
     positions = IO_stringPos(line,maxNchunks)
     tag = IO_lc(IO_stringValue(line,positions,1_pInt))        ! extract key
     select case(tag)
       case ('elasticity')
         phase_elasticity(section) = IO_lc(IO_stringValue(line,positions,2_pInt))
         do s = 1_pInt,section
           if (phase_elasticity(s) == phase_elasticity(section)) &
             phase_elasticityInstance(section) = phase_elasticityInstance(section) + 1_pInt  ! count instances
         enddo
       case ('plasticity')
         phase_plasticity(section) = IO_lc(IO_stringValue(line,positions,2_pInt))
         do s = 1_pInt,section
           if (phase_plasticity(s) == phase_plasticity(section)) &
             phase_plasticityInstance(section) = phase_plasticityInstance(section) + 1_pInt  ! count instances
         enddo
     end select
   endif
 enddo

100 end subroutine material_parsePhase


!*********************************************************************
subroutine material_parseTexture(myFile,myPart)
!*********************************************************************

 use IO
 use math, only: inRad, math_sampleRandomOri
 
 implicit none
 character(len=*), intent(in) :: myPart
 integer(pInt),    intent(in) :: myFile
 
 integer(pInt), parameter     :: maxNchunks = 13_pInt
 
 integer(pInt), dimension(1+2*maxNchunks) :: positions
 integer(pInt) :: Nsections, section, gauss, fiber, i
 character(len=64) :: tag
 character(len=1024) :: line
 
 
 Nsections = IO_countSections(myFile,myPart)
 material_Ntexture = Nsections
 if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart)

 allocate(texture_name(Nsections));     texture_name = ''
 allocate(texture_ODFfile(Nsections));  texture_ODFfile = ''
 allocate(texture_symmetry(Nsections)); texture_symmetry = 1_pInt
 allocate(texture_Ngauss(Nsections));   texture_Ngauss = 0_pInt
 allocate(texture_Nfiber(Nsections));   texture_Nfiber = 0_pInt

 texture_Ngauss = IO_countTagInPart(myFile,myPart,'(gauss)', Nsections) + &
                  IO_countTagInPart(myFile,myPart,'(random)',Nsections)
 texture_Nfiber = IO_countTagInPart(myFile,myPart,'(fiber)', Nsections)
 texture_maxNgauss = maxval(texture_Ngauss)
 texture_maxNfiber = maxval(texture_Nfiber)
 allocate(texture_Gauss   (5,texture_maxNgauss,Nsections)); texture_Gauss    = 0.0_pReal
 allocate(texture_Fiber   (6,texture_maxNfiber,Nsections)); texture_Fiber    = 0.0_pReal
 
 rewind(myFile)
 line = ''
 section = 0_pInt
 
 do while (IO_lc(IO_getTag(line,'<','>')) /= myPart)      ! wind forward to myPart
   read(myFile,'(a1024)',END=100) line
 enddo

 do
   read(myFile,'(a1024)',END=100) line
   if (IO_isBlank(line)) cycle                            ! skip empty lines
   if (IO_getTag(line,'<','>') /= '') exit                ! stop at next part
   if (IO_getTag(line,'[',']') /= '') then                ! next section
     section = section + 1_pInt
     gauss = 0_pInt
     fiber = 0_pInt
     texture_name(section) = IO_getTag(line,'[',']')
   endif
   if (section > 0_pInt) then
     positions = IO_stringPos(line,maxNchunks)
     tag = IO_lc(IO_stringValue(line,positions,1_pInt))        ! extract key
     select case(tag)

       case ('hybridia')
         texture_ODFfile(section) = IO_stringValue(line,positions,2_pInt)

       case ('symmetry')
         tag = IO_lc(IO_stringValue(line,positions,2_pInt))
         select case (tag)
           case('orthotropic')
             texture_symmetry(section) = 4_pInt
           case('monoclinic')
             texture_symmetry(section) = 2_pInt
           case default
             texture_symmetry(section) = 1_pInt
         end select
         
       case ('(random)')
         gauss = gauss + 1_pInt
         texture_Gauss(1:3,gauss,section) = math_sampleRandomOri()
         do i = 2_pInt,4_pInt,2_pInt
           tag = IO_lc(IO_stringValue(line,positions,i))
           select case (tag)
             case('scatter')
                 texture_Gauss(4,gauss,section) = IO_floatValue(line,positions,i+1_pInt)*inRad
             case('fraction')
                 texture_Gauss(5,gauss,section) = IO_floatValue(line,positions,i+1_pInt)
           end select
         enddo

       case ('(gauss)')
         gauss = gauss + 1_pInt
         do i = 2_pInt,10_pInt,2_pInt
           tag = IO_lc(IO_stringValue(line,positions,i))
           select case (tag)
             case('phi1')
                 texture_Gauss(1,gauss,section) = IO_floatValue(line,positions,i+1_pInt)*inRad
             case('phi')
                 texture_Gauss(2,gauss,section) = IO_floatValue(line,positions,i+1_pInt)*inRad
             case('phi2')
                 texture_Gauss(3,gauss,section) = IO_floatValue(line,positions,i+1_pInt)*inRad
             case('scatter')
                 texture_Gauss(4,gauss,section) = IO_floatValue(line,positions,i+1_pInt)*inRad
             case('fraction')
                 texture_Gauss(5,gauss,section) = IO_floatValue(line,positions,i+1_pInt)
           end select
         enddo

       case ('(fiber)')
         fiber = fiber + 1_pInt
         do i = 2_pInt,12_pInt,2_pInt
           tag = IO_lc(IO_stringValue(line,positions,i))
           select case (tag)
             case('alpha1')
                 texture_Fiber(1,fiber,section) = IO_floatValue(line,positions,i+1_pInt)*inRad
             case('alpha2')
                 texture_Fiber(2,fiber,section) = IO_floatValue(line,positions,i+1_pInt)*inRad
             case('beta1')
                 texture_Fiber(3,fiber,section) = IO_floatValue(line,positions,i+1_pInt)*inRad
             case('beta2')
                 texture_Fiber(4,fiber,section) = IO_floatValue(line,positions,i+1_pInt)*inRad
             case('scatter')
                 texture_Fiber(5,fiber,section) = IO_floatValue(line,positions,i+1_pInt)*inRad
             case('fraction')
                 texture_Fiber(6,fiber,section) = IO_floatValue(line,positions,i+1_pInt)
           end select
         enddo

     end select
   endif
 enddo

100 end subroutine material_parseTexture


!*********************************************************************
subroutine material_populateGrains
!*********************************************************************

 use math, only:      math_sampleRandomOri, &
                      math_sampleGaussOri, &
                      math_sampleFiberOri, &
                      math_symmetricEulers
 use mesh, only:      mesh_element, &
                      mesh_maxNips, &
                      mesh_NcpElems, &
                      mesh_ipVolume, &
                      FE_Nips
 use IO, only:        IO_error, &
                      IO_hybridIA
 use FEsolving, only: FEsolving_execIP
 use debug, only:     debug_what, &
                      debug_material, &
                      debug_levelBasic
 
 implicit none
 integer(pInt), dimension (:,:), allocatable :: Ngrains
 integer(pInt), dimension (microstructure_maxNconstituents) &
                                             :: NgrainsOfConstituent
 real(pReal), dimension (:),     allocatable :: volumeOfGrain
 real(pReal), dimension (:,:),   allocatable :: orientationOfGrain
 real(pReal), dimension (3)                  :: orientation
 real(pReal), dimension (3,3)                :: symOrientation
 integer(pInt), dimension (:),   allocatable :: phaseOfGrain, textureOfGrain
 integer(pInt) :: t,e,i,g,j,m,homog,micro,sgn,hme, myDebug
 integer(pInt) :: phaseID,textureID,dGrains,myNgrains,myNorientations, &
               grain,constituentGrain,symExtension
 real(pReal) :: extreme,rnd
 integer(pInt), dimension (:,:),   allocatable :: Nelems   ! counts number of elements in homog, micro array
 integer(pInt), dimension (:,:,:), allocatable :: elemsOfHomogMicro    ! lists element number in homog, micro array

 myDebug = debug_what(debug_material)
 
 allocate(material_volume(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ;        material_volume      = 0.0_pReal
 allocate(material_phase(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ;         material_phase       = 0_pInt
 allocate(material_texture(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ;       material_texture     = 0_pInt
 allocate(material_EulerAngles(3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; material_EulerAngles = 0.0_pReal
 
 allocate(Ngrains(material_Nhomogenization,material_Nmicrostructure)); Ngrains = 0_pInt
 allocate(Nelems(material_Nhomogenization,material_Nmicrostructure));  Nelems = 0_pInt
 
! precounting of elements for each homog/micro pair
 do e = 1_pInt, mesh_NcpElems
   homog = mesh_element(3,e)
   micro = mesh_element(4,e)
   Nelems(homog,micro) = Nelems(homog,micro) + 1_pInt
 enddo
 
 allocate(elemsOfHomogMicro(maxval(Nelems),material_Nhomogenization,material_Nmicrostructure))  
 elemsOfHomogMicro = 0_pInt

 Nelems = 0_pInt                                            ! reuse as counter

! identify maximum grain count per IP (from element) and find grains per homog/micro pair
 do e = 1_pInt,mesh_NcpElems
   homog = mesh_element(3,e)
   micro = mesh_element(4,e)
   if (homog < 1_pInt .or. homog > material_Nhomogenization) &   ! out of bounds
     call IO_error(154_pInt,e,0_pInt,0_pInt)
   if (micro < 1_pInt .or. micro > material_Nmicrostructure) &   ! out of bounds
     call IO_error(155_pInt,e,0_pInt,0_pInt)
   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
   Nelems(homog,micro)  = Nelems(homog,micro) + 1_pInt
   elemsOfHomogMicro(Nelems(homog,micro),homog,micro) = e         ! remember elements active in this homog/micro pair
   
 enddo

 allocate(volumeOfGrain(maxval(Ngrains)))           ! reserve memory for maximum case
 allocate(phaseOfGrain(maxval(Ngrains)))            ! reserve memory for maximum case
 allocate(textureOfGrain(maxval(Ngrains)))          ! reserve memory for maximum case
 allocate(orientationOfGrain(3,maxval(Ngrains)))    ! reserve memory for maximum case
 
 if (iand(myDebug,debug_levelBasic) /= 0_pInt) then
   !$OMP CRITICAL (write2out)
     write (6,*)
     write (6,*) 'MATERIAL grain population'
     write (6,*)
     write (6,'(a32,1x,a32,1x,a6)') 'homogenization_name','microstructure_name','grain#'
   !$OMP END CRITICAL (write2out)
 endif
 do homog = 1_pInt,material_Nhomogenization              ! loop over homogenizations
   dGrains = homogenization_Ngrains(homog)          ! grain number per material point
   do micro = 1_pInt,material_Nmicrostructure            ! all pairs of homog and micro
     if (Ngrains(homog,micro) > 0_pInt) then             ! an active pair of homog and micro
       myNgrains = Ngrains(homog,micro)             ! assign short name for total number of grains to populate
       if (iand(myDebug,debug_levelBasic) /= 0_pInt) then
         !$OMP CRITICAL (write2out)
           write (6,*)
           write (6,'(a32,1x,a32,1x,i6)') homogenization_name(homog),microstructure_name(micro),myNgrains
         !$OMP END CRITICAL (write2out)
       endif
     
! ----------------------------------------------------------------------------  calculate volume of each grain
       volumeOfGrain = 0.0_pReal
       grain = 0_pInt
       do hme = 1_pInt, Nelems(homog,micro)
         e = elemsOfHomogMicro(hme,homog,micro)                               ! my combination of homog and micro, only perform calculations for elements with homog, micro combinations which is indexed in cpElemsindex
         if (microstructure_elemhomo(micro)) then                             ! homogeneous distribution of grains over each element's IPs
           volumeOfGrain(grain+1_pInt:grain+dGrains) = sum(mesh_ipVolume(1:FE_Nips(mesh_element(2,e)),e))/&
                                                                         real(dGrains,pReal)
           grain = grain + dGrains                                            ! wind forward by NgrainsPerIP
         else
           forall (i = 1_pInt:FE_Nips(mesh_element(2,e))) &                        ! loop over IPs
             volumeOfGrain(grain+(i-1)*dGrains+1_pInt: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
       enddo
       
! ----------------------------------------------------------------------------  divide myNgrains as best over constituents
       NgrainsOfConstituent = 0_pInt
       forall (i = 1_pInt:microstructure_Nconstituents(micro)) &
         NgrainsOfConstituent(i) = nint(microstructure_fraction(i,micro) * myNgrains, pInt)  ! do rounding integer conversion
       do while (sum(NgrainsOfConstituent) /= myNgrains)                        ! total grain count over constituents wrong?
         sgn = sign(1_pInt, myNgrains - sum(NgrainsOfConstituent))              ! direction of required change
         extreme = 0.0_pReal
         t = 0_pInt
         do i = 1_pInt,microstructure_Nconstituents(micro)                           ! find largest deviator
           if (real(sgn,pReal)*log(NgrainsOfConstituent(i)/myNgrains/microstructure_fraction(i,micro)) > extreme) then
             extreme = real(sgn,pReal)*log(NgrainsOfConstituent(i)/myNgrains/microstructure_fraction(i,micro))
             t = i
           endif
         enddo
         NgrainsOfConstituent(t) = NgrainsOfConstituent(t) + sgn               ! change that by one
       enddo
! ----------------------------------------------------------------------------
       phaseOfGrain = 0_pInt
       textureOfGrain = 0_pInt
       orientationOfGrain = 0.0_pReal
       grain = 0_pInt                                                         ! reset microstructure grain index

       do i = 1_pInt,microstructure_Nconstituents(micro)                           ! loop over constituents
         phaseID   = microstructure_phase(i,micro)
         textureID = microstructure_texture(i,micro)
         phaseOfGrain(grain+1_pInt:grain+NgrainsOfConstituent(i)) = phaseID        ! assign resp. phase
         textureOfGrain(grain+1_pInt:grain+NgrainsOfConstituent(i)) = textureID    ! assign resp. texture

         myNorientations = ceiling(real(NgrainsOfConstituent(i),pReal)/&
                                   real(texture_symmetry(textureID),pReal),pInt)   ! max number of unique orientations (excl. symmetry)

         constituentGrain = 0_pInt                                            ! constituent grain index
                                                                              ! ---------
         if (texture_ODFfile(textureID) == '') then                           ! dealing with texture components
                                                                              ! ---------
           do t = 1_pInt,texture_Ngauss(textureID)                                 ! loop over Gauss components
             do g = 1_pInt,int(myNorientations*texture_Gauss(5,t,textureID),pInt)       ! loop over required grain count
               orientationOfGrain(:,grain+constituentGrain+g) = &
                 math_sampleGaussOri(texture_Gauss(1:3,t,textureID),&
                                     texture_Gauss(  4,t,textureID))
             enddo
             constituentGrain = constituentGrain + int(myNorientations*texture_Gauss(5,t,textureID))
           enddo

           do t = 1_pInt,texture_Nfiber(textureID)                                 ! loop over fiber components
             do g = 1_pInt,int(myNorientations*texture_Fiber(6,t,textureID),pInt)       ! loop over required grain count
               orientationOfGrain(:,grain+constituentGrain+g) = &
                 math_sampleFiberOri(texture_Fiber(1:2,t,textureID),&
                                     texture_Fiber(3:4,t,textureID),&
                                     texture_Fiber(  5,t,textureID))
             enddo
             constituentGrain = constituentGrain + int(myNorientations*texture_fiber(6,t,textureID),pInt)
           enddo

           do j = constituentGrain+1_pInt,myNorientations                          ! fill remainder with random
              orientationOfGrain(:,grain+j) = math_sampleRandomOri()
           enddo
                                                                              ! ---------
         else                                                                 ! hybrid IA
                                                                              ! ---------
           orientationOfGrain(:,grain+1:grain+myNorientations) = IO_hybridIA(myNorientations,texture_ODFfile(textureID))
           if (all(orientationOfGrain(:,grain+1) == -1.0_pReal)) call IO_error(156_pInt)  
           constituentGrain = constituentGrain + myNorientations

         endif
! ----------------------------------------------------------------------------
         symExtension = texture_symmetry(textureID) - 1_pInt
         if (symExtension > 0_pInt) then                                      ! sample symmetry
           constituentGrain = NgrainsOfConstituent(i)-myNorientations         ! calc remainder of array
           do j = 1_pInt,myNorientations                                           ! loop over each "real" orientation
             symOrientation = math_symmetricEulers(texture_symmetry(textureID),orientationOfGrain(:,j))  ! get symmetric equivalents
             e = min(symExtension,constituentGrain)                           ! are we at end of constituent grain array?
             if (e > 0_pInt) then
               orientationOfGrain(:,grain+myNorientations+1+(j-1_pInt)*symExtension:&
                                    grain+myNorientations+e+(j-1_pInt)*symExtension) = &
                 symOrientation(:,1:e)
               constituentGrain = constituentGrain - e                        ! remainder shrinks by e
             endif
           enddo
         endif
         
         grain = grain + NgrainsOfConstituent(i)                              ! advance microstructure grain index
       enddo  ! constituent

! ----------------------------------------------------------------------------
       if (.not. microstructure_elemhomo(micro)) then                           ! unless element homogeneous, reshuffle grains
         do i=1_pInt,myNgrains-1_pInt                                                     ! walk thru grains
           call random_number(rnd)
           t = nint(rnd*(myNgrains-i)+i+0.5_pReal,pInt)                         ! select a grain in remaining list
           m                       = phaseOfGrain(t)                            ! exchange current with random
           phaseOfGrain(t)         = phaseOfGrain(i)
           phaseOfGrain(i)         = m
           m                       = textureOfGrain(t)                          ! exchange current with random
           textureOfGrain(t)       = textureOfGrain(i)
           textureOfGrain(i)       = m
           orientation             = orientationOfGrain(:,t)
           orientationOfGrain(:,t) = orientationOfGrain(:,i)
           orientationOfGrain(:,i) = orientation
         enddo
       endif
       !calc fraction after weighing with volumePerGrain
       !exchange in MC steps to improve result...

! ----------------------------------------------------------------------------
       grain = 0_pInt
       do hme = 1_pInt, Nelems(homog,micro)
       e = elemsOfHomogMicro(hme,homog,micro)                                 ! only perform calculations for elements with homog, micro combinations which is indexed in cpElemsindex
         if (microstructure_elemhomo(micro)) then                             ! homogeneous distribution of grains over each element's IPs
           forall (i = 1_pInt:FE_Nips(mesh_element(2,e)), g = 1_pInt:dGrains)           ! loop over IPs and grains
             material_volume(g,i,e)        = volumeOfGrain(grain+g)
             material_phase(g,i,e)         = phaseOfGrain(grain+g)
             material_texture(g,i,e)       = textureOfGrain(grain+g)
             material_EulerAngles(:,g,i,e) = orientationOfGrain(:,grain+g)
           end forall
           FEsolving_execIP(2,e) = 1_pInt                                     ! restrict calculation to first IP only, since all other results are to be copied from this
           grain = grain + dGrains                                            ! wind forward by NgrainsPerIP
         else
           forall (i = 1_pInt:FE_Nips(mesh_element(2,e)), g = 1_pInt:dGrains)           ! loop over IPs and grains
             material_volume(g,i,e)        = volumeOfGrain(grain+(i-1_pInt)*dGrains+g)
             material_phase(g,i,e)         = phaseOfGrain(grain+(i-1_pInt)*dGrains+g)
             material_texture(g,i,e)       = textureOfGrain(grain+(i-1_pInt)*dGrains+g)
             material_EulerAngles(:,g,i,e) = orientationOfGrain(:,grain+(i-1_pInt)*dGrains+g)
           end forall
           grain = grain + FE_Nips(mesh_element(2,e)) * dGrains               ! wind forward by Nips*NgrainsPerIP
         endif
       enddo
     endif   ! active homog,micro pair
   enddo
 enddo
 
 deallocate(volumeOfGrain)
 deallocate(phaseOfGrain)
 deallocate(textureOfGrain)
 deallocate(orientationOfGrain)
 deallocate(Nelems)
 deallocate(elemsOfHomogMicro)

end subroutine material_populateGrains

end module material