From c7418db9bd39ea983c23a886924a142191e27071 Mon Sep 17 00:00:00 2001 From: Pratheek Shanthraj Date: Mon, 30 Mar 2015 09:45:10 +0000 Subject: [PATCH] FEM : hierarchical ordering of h5 output data and more meaningful visualisation for multiple grains/crystallites/phases --- code/CPFEM.f90 | 2 + code/homogenization.f90 | 120 +++++++++++++++++++++++++++++++++++++++- code/prec.f90 | 10 ++++ 3 files changed, 131 insertions(+), 1 deletion(-) diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index bbe70c1f5..9368af81f 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -343,7 +343,9 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature, dt, elFE, ip) materialpoint_F0, & materialpoint_P, & materialpoint_dPdF, & +#if defined(Marc4DAMASK) || defined(Abaqus) materialpoint_results, & +#endif materialpoint_sizeResults, & materialpoint_stressAndItsTangent, & materialpoint_postResults, & diff --git a/code/homogenization.f90 b/code/homogenization.f90 index 1b910a240..38e949b79 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -8,6 +8,9 @@ !-------------------------------------------------------------------------------------------------- module homogenization use prec, only: & +#ifdef FEM + tOutputData, & +#endif pInt, & pReal @@ -21,8 +24,16 @@ module homogenization materialpoint_P !< first P--K stress of IP real(pReal), dimension(:,:,:,:,:,:), allocatable, public :: & materialpoint_dPdF !< tangent of first P--K stress at IP +#ifdef FEM + type(tOutputData), dimension(:), allocatable, public :: & + homogOutput + type(tOutputData), dimension(:,:), allocatable, public :: & + crystalliteOutput, & + phaseOutput +#else real(pReal), dimension(:,:,:), allocatable, public :: & materialpoint_results !< results array of material point +#endif integer(pInt), public, protected :: & materialpoint_sizeResults, & homogenization_maxSizePostResults, & @@ -48,7 +59,7 @@ module homogenization damage_ID, & vacancy_concentration_ID end enum - integer(pInt), dimension(:), allocatable, private, protected :: & + integer(pInt), dimension(:), allocatable, public, protected :: & field_sizePostResults integer(pInt), dimension(:,:), allocatable, private :: & field_sizePostResult @@ -123,6 +134,9 @@ subroutine homogenization_init() constitutive_thermal_maxSizePostResults, & constitutive_vacancy_maxSizePostResults use crystallite, only: & +#ifdef FEM + crystallite_sizePostResults, & +#endif crystallite_maxSizePostResults use material use homogenization_none @@ -307,6 +321,31 @@ subroutine homogenization_init() homogenization_maxSizePostResults = max(homogenization_maxSizePostResults,homogState(p)%sizePostResults) field_maxSizePostResults = max(field_maxSizePostResults,field_sizePostResults(p)) enddo + +#ifdef FEM + allocate(homogOutput (material_Nhomogenization )) + allocate(crystalliteOutput(material_Ncrystallite, homogenization_maxNgrains)) + allocate(phaseOutput (material_Nphase, homogenization_maxNgrains)) + do p = 1, material_Nhomogenization + homogOutput(p)%sizeResults = homogState(p)%sizePostResults + field_sizePostResults(p) + homogOutput(p)%sizeIpCells = count(material_homog==p) + allocate(homogOutput(p)%output(homogOutput(p)%sizeResults,homogOutput(p)%sizeIpCells)) + enddo + do p = 1, material_Ncrystallite; do e = 1, homogenization_maxNgrains + crystalliteOutput(p,e)%sizeResults = crystallite_sizePostResults(p) + crystalliteOutput(p,e)%sizeIpCells = count(microstructure_crystallite(mesh_element(4,:)) == p .and. & + homogenization_Ngrains (mesh_element(3,:)) >= e)*mesh_maxNips + allocate(crystalliteOutput(p,e)%output(crystalliteOutput(p,e)%sizeResults,crystalliteOutput(p,e)%sizeIpCells)) + enddo; enddo + do p = 1, material_Nphase; do e = 1, homogenization_maxNgrains + phaseOutput(p,e)%sizeResults = plasticState(p)%sizePostResults + & + damageState (p)%sizePostResults + & + thermalState(p)%sizePostResults + & + vacancyState(p)%sizePostResults + phaseOutput(p,e)%sizeIpCells = count(material_phase(e,:,:) == p) + allocate(phaseOutput(p,e)%output(phaseOutput(p,e)%sizeResults,phaseOutput(p,e)%sizeIpCells)) + enddo; enddo +#else materialpoint_sizeResults = 1 & ! grain count + 1 + homogenization_maxSizePostResults & ! homogSize & homogResult + field_maxSizePostResults & ! field size & field result @@ -316,6 +355,7 @@ subroutine homogenization_init() + constitutive_thermal_maxSizePostResults & + constitutive_vacancy_maxSizePostResults) allocate(materialpoint_results(materialpoint_sizeResults,mesh_maxNips,mesh_NcpElems)) +#endif mainProcess: if (worldrank == 0) then write(6,'(/,a)') ' <<<+- homogenization init -+>>>' @@ -342,7 +382,9 @@ subroutine homogenization_init() write(6,'(a32,1x,7(i8,1x))') 'materialpoint_requested: ', shape(materialpoint_requested) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_converged: ', shape(materialpoint_converged) write(6,'(a32,1x,7(i8,1x),/)') 'materialpoint_doneAndHappy: ', shape(materialpoint_doneAndHappy) +#ifndef FEM write(6,'(a32,1x,7(i8,1x),/)') 'materialpoint_results: ', shape(materialpoint_results) +#endif write(6,'(a32,1x,7(i8,1x))') 'maxSizePostResults: ', homogenization_maxSizePostResults endif flush(6) @@ -703,9 +745,18 @@ subroutine materialpoint_postResults FEsolving_execElem, & FEsolving_execIP use mesh, only: & +#ifdef FEM + mesh_maxNips, & +#endif mesh_element use material, only: & mappingHomogenization, & +#ifdef FEM + mappingConstitutive, & + homogenization_maxNgrains, & + material_Ncrystallite, & + material_Nphase, & +#endif homogState, & plasticState, & damageState, & @@ -715,8 +766,17 @@ subroutine materialpoint_postResults homogenization_Ngrains, & microstructure_crystallite use constitutive, only: & +#ifdef FEM + constitutive_maxSizePostResults, & + constitutive_damage_maxSizePostResults, & + constitutive_thermal_maxSizePostResults, & + constitutive_vacancy_maxSizePostResults, & +#endif constitutive_postResults use crystallite, only: & +#ifdef FEM + crystallite_maxSizePostResults, & +#endif crystallite_sizePostResults, & crystallite_postResults @@ -729,6 +789,63 @@ subroutine materialpoint_postResults g, & !< grain number i, & !< integration point number e !< element number +#ifdef FEM + integer(pInt) :: & + myHomog, & + myPhase, & + crystalliteCtr(material_Ncrystallite, homogenization_maxNgrains), & + phaseCtr (material_Nphase, homogenization_maxNgrains) + real(pReal), dimension(1+crystallite_maxSizePostResults + & + 1+constitutive_maxSizePostResults + & + constitutive_damage_maxSizePostResults + & + constitutive_thermal_maxSizePostResults + & + constitutive_vacancy_maxSizePostResults) :: & + crystalliteResults + + + + crystalliteCtr = 0_pInt; phaseCtr = 0_pInt + elementLooping: do e = FEsolving_execElem(1),FEsolving_execElem(2) + myNgrains = homogenization_Ngrains(mesh_element(3,e)) + myCrystallite = microstructure_crystallite(mesh_element(4,e)) + IpLooping: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) + myHomog = mappingHomogenization(2,i,e) + thePos = mappingHomogenization(1,i,e) + homogOutput(myHomog)%output(1: & + homogState(myHomog)%sizePostResults, & + thePos) = homogenization_postResults(i,e) + homogOutput(myHomog)%output(homogState(myHomog)%sizePostResults+1: & + homogState(myHomog)%sizePostResults+field_sizePostResults(myHomog), & + thePos) = field_postResults(i,e) + + grainLooping :do g = 1,myNgrains + myPhase = mappingConstitutive(2,g,i,e) + crystalliteResults(1:1+crystallite_sizePostResults(myCrystallite) + & + 1+plasticState(myPhase)%sizePostResults + & + damageState (myPhase)%sizePostResults + & + thermalState(myPhase)%sizePostResults + & + vacancyState(myPhase)%sizePostResults) = crystallite_postResults(g,i,e) + if (microstructure_crystallite(mesh_element(4,e)) == myCrystallite .and. & + homogenization_Ngrains (mesh_element(3,e)) >= g) then + crystalliteCtr(myCrystallite,g) = crystalliteCtr(myCrystallite,g) + 1_pInt + crystalliteOutput(myCrystallite,g)% & + output(1:crystalliteOutput(myCrystallite,g)%sizeResults,crystalliteCtr(myCrystallite,g)) = & + crystalliteResults(2:1+crystallite_sizePostResults(myCrystallite)) + endif + if (material_phase(g,i,e) == myPhase) then + phaseCtr(myPhase,g) = phaseCtr(myPhase,g) + 1_pInt + phaseOutput(myPhase,g)% & + output(1:phaseOutput(myPhase,g)%sizeResults,phaseCtr(myPhase,g)) = & + crystalliteResults(3 + crystallite_sizePostResults(myCrystallite): & + plasticState(myphase)%sizePostResults + & + damageState (myphase)%sizePostResults + & + thermalState(myphase)%sizePostResults + & + vacancyState(myphase)%sizePostResults) + endif + enddo grainLooping + enddo IpLooping + enddo elementLooping +#else !$OMP PARALLEL DO PRIVATE(myNgrains,myCrystallite,thePos,theSize) elementLooping: do e = FEsolving_execElem(1),FEsolving_execElem(2) @@ -767,6 +884,7 @@ subroutine materialpoint_postResults enddo IpLooping enddo elementLooping !$OMP END PARALLEL DO +#endif end subroutine materialpoint_postResults diff --git a/code/prec.f90 b/code/prec.f90 index 3e69b9350..e86b6f277 100644 --- a/code/prec.f90 +++ b/code/prec.f90 @@ -101,6 +101,16 @@ type, public :: p_intvec field !< field data end type +#ifdef FEM + type, public :: tOutputData + integer(pInt) :: & + sizeIpCells = 0_pInt , & + sizeResults = 0_pInt + real(pReal), allocatable, dimension(:,:) :: & + output !< output data + end type +#endif + public :: & prec_init