From 84c43741a6d4f8a9e188645c0b475897ab4466c6 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 29 May 2013 17:23:49 +0000 Subject: [PATCH] reworked distribution of microstructure constituents. now each material point represents "as good as possible" the fractional content of constituents. removed error for volume fractions not equalling 1. implemented capability to rotate the texture given in material.config. --- code/material.f90 | 298 +++++++++++++++++++++++++++++----------------- 1 file changed, 189 insertions(+), 109 deletions(-) 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)