Calculation of current ip volumes now working. Crystallite output also reflects current grain volume, not reference volume. However, this is only available for Marc. Abaqus and spectral method still return the reference ip volume. The ip coordinates though are correctly updated for all solver types.

This commit is contained in:
Christoph Kords 2011-08-10 16:37:17 +00:00
parent 34de2e301b
commit 1ffb59a96a
3 changed files with 55 additions and 32 deletions

View File

@ -275,6 +275,7 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
mesh_ipCenterOfGravity, &
mesh_build_subNodeCoords, &
mesh_build_ipVolumes, &
mesh_build_ipCoordinates, &
FE_Nips, &
FE_Nnodes
use material, only: homogenization_maxNgrains, &
@ -532,9 +533,10 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
write(6,'(a,i8,a,i8)') '<< CPFEM >> Calculation for elements ',FEsolving_execElem(1),' to ',FEsolving_execElem(2)
!$OMP END CRITICAL (write2out)
endif
if (any(.not. crystallite_localConstitution) .and. FEsolver == 'Marc') then
if (FEsolver == 'Marc') then ! marc updates nodal coordinates, whereas Abaqus and spectral solver directly update ip coordinates. In the latter case it is not possible to get the current ip volume, since the current nodal positions are unknown
call mesh_build_subNodeCoords() ! update subnodal coordinates
call mesh_build_ipVolumes() ! update ip center of gravity
call mesh_build_ipCoordinates() ! update ip coordinates
call mesh_build_ipVolumes() ! update ip volumes
endif
if (debug_verbosity > 0) then
!$OMP CRITICAL (write2out)

View File

@ -3066,12 +3066,13 @@ function crystallite_postResults(&
math_I3, &
inDeg, &
math_Mandel6to33
use mesh, only: mesh_element
use mesh, only: mesh_element, &
mesh_ipVolume
use material, only: microstructure_crystallite, &
crystallite_Noutput, &
material_phase, &
material_texture, &
material_volume
homogenization_Ngrains
use constitutive, only: constitutive_sizePostResults, &
constitutive_postResults
@ -3109,7 +3110,7 @@ function crystallite_postResults(&
crystallite_postResults(c+1) = material_texture(g,i,e) ! textureID of grain
case ('volume')
mySize = 1_pInt
crystallite_postResults(c+1) = material_volume(g,i,e) ! grain volume (not fraction but absolute, right?)
crystallite_postResults(c+1) = mesh_ipVolume(i,e) / homogenization_Ngrains(mesh_element(3,e)) ! grain volume (not fraction but absolute)
case ('orientation')
mySize = 4_pInt
crystallite_postResults(c+1:c+mySize) = crystallite_orientation(1:4,g,i,e) ! grain orientation as quaternion

View File

@ -322,6 +322,7 @@
close (fileUnit)
call mesh_build_subNodeCoords()
call mesh_build_ipCoordinates()
call mesh_build_ipVolumes()
call mesh_build_ipAreas()
call mesh_build_nodeTwins()
@ -3134,32 +3135,23 @@ endsubroutine
!***********************************************************
! calculation of IP volume
! calculation of IP coordinates
!
! allocate globals
! _ipVolume
! _ipCenterOfGravity
!***********************************************************
subroutine mesh_build_ipVolumes()
subroutine mesh_build_ipCoordinates()
use prec, only: pInt, tol_gravityNodePos
use math, only: math_volTetrahedron
implicit none
integer(pInt) e,f,t,i,j,k,n
integer(pInt), parameter :: Ntriangles = FE_NipFaceNodes-2 ! each interface is made up of this many triangles
logical(pInt), dimension(mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNode ! flagList to find subnodes determining center of grav
real(pReal), dimension(3,mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNodePos ! coordinates of subnodes determining center of grav
real(pReal), dimension(3,FE_NipFaceNodes) :: nPos ! coordinates of nodes on IP face
real(pReal), dimension(Ntriangles,FE_NipFaceNodes) :: volume ! volumes of possible tetrahedra
real(pReal), dimension(3) :: centerOfGravity
logical :: calcIPvolume = .false.
if (.not. allocated(mesh_ipVolume)) then
allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems))
if (.not. allocated(mesh_ipCenterOfGravity)) then
allocate(mesh_ipCenterOfGravity(3,mesh_maxNips,mesh_NcpElems))
mesh_ipVolume = 0.0_pReal
mesh_ipCenterOfGravity = 0.0_pReal
calcIPvolume = .true.
endif
do e = 1,mesh_NcpElems ! loop over cpElems
@ -3187,23 +3179,51 @@ endsubroutine
enddo
centerOfGravity = sum(gravityNodePos,2)/count(gravityNode)
mesh_ipCenterOfGravity(:,i,e) = centerOfGravity
if (calcIPvolume) then
do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP and add tetrahedra which connect to CoG
forall (n = 1:FE_NipFaceNodes) nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e)
forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles) & ! start at each interface node and build valid triangles to cover interface
volume(j,n) = math_volTetrahedron(nPos(:,n), & ! calc volume of respective tetrahedron to CoG
nPos(:,1+mod(n-1 +j ,FE_NipFaceNodes)), & ! start at offset j
nPos(:,1+mod(n-1 +j+1,FE_NipFaceNodes)), & ! and take j's neighbor
centerOfGravity)
mesh_ipVolume(i,e) = mesh_ipVolume(i,e) + sum(volume) ! add contribution from this interface
enddo
mesh_ipVolume(i,e) = mesh_ipVolume(i,e) / FE_NipFaceNodes ! renormalize with interfaceNodeNum due to loop over them
endif
enddo
enddo
endsubroutine
!***********************************************************
! calculation of IP volume
!
! allocate globals
! _ipVolume
!***********************************************************
subroutine mesh_build_ipVolumes()
use prec, only: pInt, tol_gravityNodePos
use math, only: math_volTetrahedron
implicit none
integer(pInt) e,f,t,i,j,n
integer(pInt), parameter :: Ntriangles = FE_NipFaceNodes-2 ! each interface is made up of this many triangles
real(pReal), dimension(3,FE_NipFaceNodes) :: nPos ! coordinates of nodes on IP face
real(pReal), dimension(Ntriangles,FE_NipFaceNodes) :: volume ! volumes of possible tetrahedra
if (.not. allocated(mesh_ipVolume)) then
allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems))
endif
mesh_ipVolume = 0.0_pReal
do e = 1,mesh_NcpElems ! loop over cpElems
t = mesh_element(2,e) ! get elemType
do i = 1,FE_Nips(t) ! loop over IPs of elem
do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP and add tetrahedra which connect to CoG
forall (n = 1:FE_NipFaceNodes) &
nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e)
forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles) & ! start at each interface node and build valid triangles to cover interface
volume(j,n) = math_volTetrahedron(nPos(:,n), & ! calc volume of respective tetrahedron to CoG
nPos(:,1+mod(n-1 +j ,FE_NipFaceNodes)), & ! start at offset j
nPos(:,1+mod(n-1 +j+1,FE_NipFaceNodes)), & ! and take j's neighbor
mesh_ipCenterOfGravity(:,i,e))
mesh_ipVolume(i,e) = mesh_ipVolume(i,e) + sum(volume) ! add contribution from this interface
enddo
mesh_ipVolume(i,e) = mesh_ipVolume(i,e) / FE_NipFaceNodes ! renormalize with interfaceNodeNum due to loop over them
enddo
enddo
endsubroutine