diff --git a/code/material.f90 b/code/material.f90 index 5b91a5de1..dfb9c0178 100644 --- a/code/material.f90 +++ b/code/material.f90 @@ -116,7 +116,7 @@ module material material_volume, & !< volume of each grain,IP,element texture_Gauss, & !< data of each Gauss component texture_Fiber, & !< data of each Fiber component - texture_Rotation !< rotation of each texture + texture_rotation !< rotation of each texture logical, dimension(:), allocatable, private :: & homogenization_active @@ -196,22 +196,22 @@ subroutine material_init if (minval(microstructure_texture(1:microstructure_Nconstituents(m),m)) < 1_pInt .or. & maxval(microstructure_texture(1:microstructure_Nconstituents(m),m)) > material_Ntexture) & call IO_error(152_pInt,m) - if (abs(sum(microstructure_fraction(:,m)) - 1.0_pReal) >= 1.0e-10_pReal) then - if (iand(myDebug,debug_levelExtensive) /= 0_pInt) then - write(6,'(a)') ' sum of microstructure fraction = ',sum(microstructure_fraction(:,m)) - endif - call IO_error(153_pInt,m) - endif +! if (abs(sum(microstructure_fraction(:,m)) - 1.0_pReal) >= 1.0e-6_pReal) then ! have ppm precision in fractions +! if (iand(myDebug,debug_levelExtensive) /= 0_pInt) then +! write(6,'(a,1x,f12.9)') ' sum of microstructure fraction = ',sum(microstructure_fraction(:,m)) +! endif +! call IO_error(153_pInt,m) +! endif enddo debugOut: if (iand(myDebug,debug_levelExtensive) /= 0_pInt) then write(6,'(/,a,/)') ' MATERIAL configuration' write(6,'(a32,1x,a16,1x,a6)') 'homogenization ','type ','grains' do h = 1_pInt,material_Nhomogenization - write(6,'(1x,a32,1x,a16,1x,i4)') homogenization_name(h),homogenization_type(h),homogenization_Ngrains(h) + write(6,'(1x,a32,1x,a16,1x,i6)') homogenization_name(h),homogenization_type(h),homogenization_Ngrains(h) enddo - write(6,'(/,a32,19x,a11,1x,a12,1x,a13)') 'microstructure','crystallite','constituents','homogeneous' + write(6,'(/,a14,18x,1x,a11,1x,a12,1x,a13)') 'microstructure','crystallite','constituents','homogeneous' do m = 1_pInt,material_Nmicrostructure - write(6,'(a32,4x,i4,8x,i4,8x,l1)') microstructure_name(m), & + write(6,'(1x,a32,1x,i11,1x,i12,1x,l13)') microstructure_name(m), & microstructure_crystallite(m), & microstructure_Nconstituents(m), & microstructure_elemhomo(m) @@ -603,9 +603,9 @@ subroutine material_parseTexture(myFile,myPart) 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 - allocate(texture_Rotation(3,3,Nsections)); - do j =1_pInt, Nsections - texture_Rotation(1:3,1:3,j) = math_I3 + allocate(texture_rotation(3,3,Nsections)); + do j = 1_pInt, Nsections + texture_rotation(1:3,1:3,j) = math_I3 enddo rewind(myFile) @@ -636,23 +636,23 @@ subroutine material_parseTexture(myFile,myPart) textureType: select case(tag) case ('rotation') textureType - do j = 1_pInt, 3_pInt + do j = 1_pInt, 3_pInt ! look for "x", "y", and "z" entries tag = IO_lc(IO_stringValue(line,positions,j+1_pInt)) select case (tag) case('x', '+x') - texture_Rotation(1_pInt:3_pInt,j,section) = (/1.0_pReal, 0.0_pReal, 0.0_pReal/) ! original axis is now +x-axis + texture_rotation(j,1:3,section) = (/ 1.0_pReal, 0.0_pReal, 0.0_pReal/) ! original axis is now +x-axis case('-x') - texture_Rotation(1_pInt:3_pInt,j,section) = (/-1.0_pReal, 0.0_pReal, 0.0_pReal/) ! original axis is now -x-axis + texture_rotation(j,1:3,section) = (/-1.0_pReal, 0.0_pReal, 0.0_pReal/) ! original axis is now -x-axis case('y', '+y') - texture_Rotation(1_pInt:3_pInt,j,section) = (/0.0_pReal, 1.0_pReal, 0.0_pReal/) ! original axis is now +y-axis + texture_rotation(j,1:3,section) = (/ 0.0_pReal, 1.0_pReal, 0.0_pReal/) ! original axis is now +y-axis case('-y') - texture_Rotation(1_pInt:3_pInt,j,section) = (/0.0_pReal, -1.0_pReal, 0.0_pReal/) ! original axis is now -y-axis + texture_rotation(j,1:3,section) = (/ 0.0_pReal,-1.0_pReal, 0.0_pReal/) ! original axis is now -y-axis case('z', '+z') - texture_Rotation(1_pInt:3_pInt,j,section) = (/0.0_pReal, 0.0_pReal, 1.0_pReal/) ! original axis is now +z-axis + texture_rotation(j,1:3,section) = (/ 0.0_pReal, 0.0_pReal, 1.0_pReal/) ! original axis is now +z-axis case('-z') - texture_Rotation(1_pInt:3_pInt,j,section) = (/0.0_pReal, 0.0_pReal, -1.0_pReal/) ! original axis is now -z-axis + texture_rotation(j,1:3,section) = (/ 0.0_pReal, 0.0_pReal,-1.0_pReal/) ! original axis is now -z-axis case default - call IO_error(156_pInt) + call IO_error(157_pInt,section) end select enddo @@ -735,6 +735,10 @@ subroutine material_parseTexture(myFile,myPart) !-------------------------------------------------------------------------------------------------- subroutine material_populateGrains use math, only: & + math_RtoEuler, & + math_EulerToR, & + math_mul33x33, & + math_range, & math_sampleRandomOri, & math_sampleGaussOri, & math_sampleFiberOri, & @@ -758,16 +762,20 @@ subroutine material_populateGrains implicit none integer(pInt), dimension (:,:), allocatable :: Ngrains - integer(pInt), dimension (microstructure_maxNconstituents) & - :: NgrainsOfConstituent + integer(pInt), dimension (microstructure_maxNconstituents) :: & + NgrainsOfConstituent, & + currentGrainOfConstituent, & + randomOrder + real(pReal), dimension (microstructure_maxNconstituents) :: & + rndArray 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 + integer(pInt) :: t,e,i,g,j,m,c,r,homog,micro,sgn,hme, myDebug, & + phaseID,textureID,dGrains,myNgrains,myNorientations,myNconstituents, & + grain,constituentGrain,ipGrain,symExtension, ip real(pReal) :: extreme,rnd integer(pInt), dimension (:,:), allocatable :: Nelems ! counts number of elements in homog, micro array type(p_intvec), dimension (:,:), allocatable :: elemsOfHomogMicro ! lists element number in homog, micro array @@ -803,23 +811,23 @@ subroutine material_populateGrains ! identify maximum grain count per IP (from element) and find grains per homog/micro pair Nelems = 0_pInt ! reuse as counter elementLooping: do e = 1_pInt,mesh_NcpElems - t = FE_geomtype(mesh_element(2,e)) + t = FE_geomtype(mesh_element(2,e)) 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) + if (microstructure_elemhomo(micro)) then ! how many grains are needed at this element? + dGrains = homogenization_Ngrains(homog) ! only one set of Ngrains (other IPs are plain copies) else - dGrains = homogenization_Ngrains(homog) * FE_Nips(t) + dGrains = homogenization_Ngrains(homog) * FE_Nips(t) ! each IP has Ngrains endif - Ngrains(homog,micro) = Ngrains(homog,micro) + dGrains - Nelems(homog,micro) = Nelems(homog,micro) + 1_pInt + Ngrains(homog,micro) = Ngrains(homog,micro) + dGrains ! total grain count + Nelems(homog,micro) = Nelems(homog,micro) + 1_pInt ! total element count elemsOfHomogMicro(homog,micro)%p(Nelems(homog,micro)) = e ! remember elements active in this homog/micro pair - enddo elementLooping + 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 @@ -836,41 +844,58 @@ subroutine material_populateGrains 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 + myNconstituents = microstructure_Nconstituents(micro) ! assign short name for number of constituents if (iand(myDebug,debug_levelBasic) /= 0_pInt) then !$OMP CRITICAL (write2out) 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(homog,micro)%p(hme) ! my combination of homog and micro, only perform calculations for elements with homog, micro combinations which is indexed in cpElemsindex + e = elemsOfHomogMicro(homog,micro)%p(hme) ! my combination of homog and micro, only perform calculations for elements with homog, micro combinations which is indexed in cpElemsindex t = FE_geomtype(mesh_element(2,e)) 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(t),e))/& - real(dGrains,pReal) - grain = grain + dGrains ! wind forward by NgrainsPerIP + real(dGrains,pReal) ! each grain combines size of all IPs in that element + grain = grain + dGrains ! wind forward by Ngrains@IP else forall (i = 1_pInt:FE_Nips(t)) & ! 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(t) * dGrains ! wind forward by Nips*NgrainsPerIP + mesh_ipVolume(i,e)/dGrains ! assign IPvolume/Ngrains@IP to all grains of IP + grain = grain + FE_Nips(t) * dGrains ! wind forward by Nips*Ngrains@IP endif enddo - + + if (grain /= myNgrains) & + call IO_error(0,e = homog,i = micro,ext_msg = 'inconsistent grain count after volume calc') + !-------------------------------------------------------------------------------------------------- ! divide myNgrains as best over constituents - NgrainsOfConstituent = 0_pInt - forall (i = 1_pInt:microstructure_Nconstituents(micro)) & +! +! example: three constituents with fractions of 0.25, 0.25, and 0.5 distributed over 20 (microstructure) grains +! +! ***** ***** ********** +! NgrainsOfConstituent: 5, 5, 10 +! counters: +! |-----> grain (if constituent == 2) +! |--> constituentGrain (of constituent 2) +! + + NgrainsOfConstituent = 0_pInt ! reset counter of grains per constituent + forall (i = 1_pInt:myNconstituents) & 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 + do i = 1_pInt,myNconstituents ! 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 @@ -879,116 +904,171 @@ subroutine material_populateGrains NgrainsOfConstituent(t) = NgrainsOfConstituent(t) + sgn ! change that by one enddo +!-------------------------------------------------------------------------------------------------- +! assign phase and texture info + phaseOfGrain = 0_pInt textureOfGrain = 0_pInt orientationOfGrain = 0.0_pReal - grain = 0_pInt ! reset microstructure grain index - constituents: do i = 1_pInt,microstructure_Nconstituents(micro) ! loop over constituents + texture: do i = 1_pInt,myNconstituents ! loop over constituents + grain = sum(NgrainsOfConstituent(1_pInt:i-1_pInt)) ! set microstructure grain index of current constituent + ! "grain" points to start of this constituent's grain population + constituentGrain = 0_pInt ! constituent grain index + phaseID = microstructure_phase(i,micro) textureID = microstructure_texture(i,micro) - phaseOfGrain(grain+1_pInt:grain+NgrainsOfConstituent(i)) = phaseID ! assign resp. phase + 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 - !-------------------------------------------------------------------------------------------------- -! dealing with texture components +! ...has texture components if (texture_ODFfile(textureID) == '') then - do t = 1_pInt,texture_Ngauss(textureID) ! loop over Gauss components + gauss: 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 + constituentGrain = & + constituentGrain + int(myNorientations*texture_Gauss(5,t,textureID)) ! advance counter for grains of current constituent + enddo gauss - do t = 1_pInt,texture_Nfiber(textureID) ! loop over fiber components + fiber: 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 + constituentGrain = & + constituentGrain + int(myNorientations*texture_fiber(6,t,textureID),pInt) ! advance counter for grains of current constituent + enddo fiber + random: do constituentGrain = constituentGrain+1_pInt,myNorientations ! fill remainder with random + orientationOfGrain(:,grain+constituentGrain) = math_sampleRandomOri() + enddo random !-------------------------------------------------------------------------------------------------- -! 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 +! ...has hybrid IA + else + orientationOfGrain(1:3,grain+1_pInt:grain+myNorientations) = & + IO_hybridIA(myNorientations,texture_ODFfile(textureID)) + if (all(orientationOfGrain(1:3,grain+1_pInt) == -1.0_pReal)) call IO_error(156_pInt) endif +!-------------------------------------------------------------------------------------------------- +! ...texture rotation + + do j = 1_pInt,myNorientations ! loop over each "real" orientation + orientationOfGrain(1:3,grain+j) = math_RtoEuler( & ! translate back to Euler angles + math_mul33x33( & ! pre-multiply + math_EulertoR(orientationOfGrain(1:3,grain+j)), & ! face-value orientation + texture_rotation(1:3,1:3,textureID) & ! rotation matrix and + ) & + ) + enddo + +!-------------------------------------------------------------------------------------------------- +! ...sample symmetry + symExtension = texture_symmetry(textureID) - 1_pInt - if (symExtension > 0_pInt) then ! sample symmetry - constituentGrain = NgrainsOfConstituent(i)-myNorientations ! calc remainder of array + if (symExtension > 0_pInt) then ! sample symmetry (number of additional equivalent orientations) + constituentGrain = myNorientations ! start right after "real" orientations 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? + symOrientation = math_symmetricEulers(texture_symmetry(textureID), & + orientationOfGrain(1:3,grain+j)) ! get symmetric equivalents + e = min(symExtension,NgrainsOfConstituent(i)-constituentGrain) ! do not overshoot 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 + orientationOfGrain(1:3,grain+constituentGrain+1: & + grain+constituentGrain+e) = & + symOrientation(1:3,1:e) + constituentGrain = constituentGrain + e ! remainder shrinks by e endif enddo endif - - grain = grain + NgrainsOfConstituent(i) ! advance microstructure grain index - enddo constituents !-------------------------------------------------------------------------------------------------- -! unless element homogeneous, reshuffle grains - if (.not. microstructure_elemhomo(micro)) then - do i=1_pInt,myNgrains-1_pInt ! walk thru grains +! shuffle grains within current constituent + + do j = 1_pInt,NgrainsOfConstituent(i)-1_pInt ! walk thru grains of current constituent 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 + t = nint(rnd*(NgrainsOfConstituent(i)-j)+j+0.5_pReal,pInt) ! select a grain in remaining list + m = phaseOfGrain(grain+t) ! exchange current with random + phaseOfGrain(grain+t) = phaseOfGrain(grain+j) + phaseOfGrain(grain+j) = m + m = textureOfGrain(grain+t) ! exchange current with random + textureOfGrain(grain+t) = textureOfGrain(grain+j) + textureOfGrain(grain+j) = m + orientation = orientationOfGrain(1:3,grain+t) ! exchange current with random + orientationOfGrain(1:3,grain+t) = orientationOfGrain(1:3,grain+j) + orientationOfGrain(1:3,grain+j) = orientation enddo - endif - + + enddo texture + !-------------------------------------------------------------------------------------------------- -! calc fraction after weighing with volumePerGrain, exchange in MC steps to improve result - grain = 0_pInt +!>@ToDo calc fraction after weighing with volumePerGrain, exchange in MC steps to improve result (humbug at the moment) + +!-------------------------------------------------------------------------------------------------- +! distribute grains of all constituents as accurately as possible to given constituent fractions + + ip = 0_pInt + currentGrainOfConstituent = 0_pInt + do hme = 1_pInt, Nelems(homog,micro) e = elemsOfHomogMicro(homog,micro)%p(hme) ! only perform calculations for elements with homog, micro combinations which is indexed in cpElemsindex t = FE_geomtype(mesh_element(2,e)) if (microstructure_elemhomo(micro)) then ! homogeneous distribution of grains over each element's IPs - forall (i = 1_pInt:FE_Nips(t), 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 + m = 1_pInt ! process only first IP else - forall (i = 1_pInt:FE_Nips(t), 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(t) * dGrains ! wind forward by Nips*NgrainsPerIP + m = FE_Nips(t) ! process all IPs endif + + do i = 1_pInt, m ! loop over necessary IPs + ip = ip + 1_pInt ! keep track of total ip count + ipGrain = 0_pInt ! count number of grains assigned at this IP + randomOrder = math_range(microstructure_maxNconstituents) ! start out with ordered sequence of constituents + call random_number(rndArray) ! as many rnd numbers as (max) constituents + do j = 1_pInt, myNconstituents - 1_pInt ! loop over constituents ... + r = nint(rndArray(j)*(myNconstituents-j)+j+0.5_pReal,pInt) ! ... select one in remaining list + c = randomOrder(r) ! ... call it "c" + randomOrder(r) = randomOrder(j) ! ... and exchange with present position in constituent list + grain = sum(NgrainsOfConstituent(1:c-1_pInt)) ! figure out actual starting index in overall/consecutive grain population + do g = 1_pInt, min(dGrains-ipGrain, & ! leftover number of grains at this IP + max(0_pInt, & ! no negative values + nint(real(ip * dGrains * NgrainsOfConstituent(c)) / & ! fraction of grains scaled to this constituent... + real(myNgrains),pInt) - & ! ...minus those already distributed + currentGrainOfConstituent(c))) + ipGrain = ipGrain + 1_pInt ! advance IP grain counter + currentGrainOfConstituent(c) = currentGrainOfConstituent(c) + 1_pInt ! advance index of grain population for constituent c + material_volume(ipGrain,i,e) = volumeOfGrain(grain+currentGrainOfConstituent(c)) ! assign properties + material_phase(ipGrain,i,e) = phaseOfGrain(grain+currentGrainOfConstituent(c)) + material_texture(ipGrain,i,e) = textureOfGrain(grain+currentGrainOfConstituent(c)) + material_EulerAngles(1:3,ipGrain,i,e) = orientationOfGrain(1:3,grain+currentGrainOfConstituent(c)) + enddo; enddo + + c = randomOrder(microstructure_Nconstituents(micro)) ! look up constituent remaining after random shuffling + grain = sum(NgrainsOfConstituent(1:c-1_pInt)) ! figure out actual starting index in overall/consecutive grain population + do ipGrain = ipGrain + 1_pInt, dGrains ! ensure last constituent fills up to dGrains + currentGrainOfConstituent(c) = currentGrainOfConstituent(c) + 1_pInt + material_volume(ipGrain,i,e) = volumeOfGrain(grain+currentGrainOfConstituent(c)) + material_phase(ipGrain,i,e) = phaseOfGrain(grain+currentGrainOfConstituent(c)) + material_texture(ipGrain,i,e) = textureOfGrain(grain+currentGrainOfConstituent(c)) + material_EulerAngles(1:3,ipGrain,i,e) = orientationOfGrain(1:3,grain+currentGrainOfConstituent(c)) + enddo + + enddo + + do i = i, FE_Nips(t) ! loop over IPs to (possibly) distribute copies from first IP + material_volume (1_pInt:dGrains,i,e) = material_volume (1_pInt:dGrains,1,e) + material_phase (1_pInt:dGrains,i,e) = material_phase (1_pInt:dGrains,1,e) + material_texture(1_pInt:dGrains,i,e) = material_texture(1_pInt:dGrains,1,e) + material_EulerAngles(1:3,1_pInt:dGrains,i,e) = material_EulerAngles(1:3,1_pInt:dGrains,1,e) + enddo + enddo endif ! active homog,micro pair enddo @@ -999,10 +1079,10 @@ subroutine material_populateGrains deallocate(textureOfGrain) deallocate(orientationOfGrain) deallocate(Nelems) - !> ToDo - causing segmentation fault: needs looking into + !>@ToDo - causing segmentation fault: needs looking into !do homog = 1,material_Nhomogenization ! do micro = 1,material_Nmicrostructure - ! if (Nelems(homog,micro) > 0_pInt) deallocate(elemsOfHomogMicro(homog,micro)%p) + ! if (Nelems(homog,micro) > 0_pInt) deallocate(elemsOfHomogMicro(homog,micro)%p) ! enddo !enddo deallocate(elemsOfHomogMicro)