Merge branch 'marc-mesh-simplification' of magit1.mpie.de:/damask/DAMASK into marc-mesh-simplification

This commit is contained in:
Martin Diehl 2019-06-12 19:56:57 +02:00
commit 9a9e99739e
3 changed files with 187 additions and 182 deletions

View File

@ -16,6 +16,7 @@ module material
use debug
use mesh
use numerics
use rotations
use discretization
implicit none
@ -152,7 +153,7 @@ module material
damageState
integer, dimension(:,:,:), allocatable, public, protected :: &
material_texture !< texture (index) of each grain,IP,element. Used only by plastic_nonlocal
material_texture !< texture (index) of each grain,IP,element. Only used by plastic_nonlocal
real(pReal), dimension(:,:,:,:), allocatable, public, protected :: &
material_EulerAngles !< initial orientation of each grain,IP,element
@ -172,11 +173,8 @@ module material
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, private :: &
texture_Gauss, & !< data of each Gauss component
texture_transformation !< transformation for each texture
microstructure_fraction !< vol fraction of each constituent in microstructure
logical, dimension(:), allocatable, private :: &
homogenization_active
@ -561,6 +559,8 @@ subroutine material_parseMicrostructure
call IO_error(153,ext_msg=microstructure_name(m))
enddo
call config_deallocate('material.config/microstructure')
end subroutine material_parseMicrostructure
@ -707,69 +707,66 @@ end subroutine material_parsePhase
!--------------------------------------------------------------------------------------------------
subroutine material_parseTexture
integer :: section, gauss, j, t, i
character(len=65536), dimension(:), allocatable :: strings ! Values for given key in material config
integer, dimension(:), allocatable :: chunkPos
integer :: j, t, i
character(len=65536), dimension(:), allocatable :: strings ! Values for given key in material config
integer, dimension(:), allocatable :: chunkPos
real(pReal), dimension(3,3) :: texture_transformation ! maps texture to microstructure coordinate system
type(rotation) :: eulers
do t=1, size(config_texture)
if (config_texture(t)%countKeys('(gauss)') /= 1) call IO_error(147,ext_msg='count((gauss)) != 1')
if (config_texture(t)%keyExists('symmetry')) call IO_error(147,ext_msg='symmetry')
if (config_texture(t)%keyExists('(random)')) call IO_error(147,ext_msg='(random)')
if (config_texture(t)%keyExists('(fiber)')) call IO_error(147,ext_msg='(fiber)')
enddo
do t=1, size(config_texture)
if (config_texture(t)%countKeys('(gauss)') /= 1) call IO_error(147,ext_msg='count((gauss)) !=1')
if (config_texture(t)%keyExists('symmetry')) call IO_error(147,ext_msg='symmetry')
if (config_texture(t)%keyExists('(random)')) call IO_error(147,ext_msg='(random)')
if (config_texture(t)%keyExists('(fiber)')) call IO_error(147,ext_msg='(fiber)')
enddo
allocate(texture_Gauss (3,size(config_texture)), source=0.0_pReal)
allocate(texture_Gauss (5,1,size(config_texture)), source=0.0_pReal)
allocate(texture_transformation(3,3,size(config_texture)), source=0.0_pReal)
texture_transformation = spread(math_I3,3,size(config_texture))
do t=1, size(config_texture)
section = t
gauss = 0
if (config_texture(t)%keyExists('axes')) then
strings = config_texture(t)%getStrings('axes')
do j = 1, 3 ! look for "x", "y", and "z" entries
select case (strings(j))
case('x', '+x')
texture_transformation(j,1:3,t) = [ 1.0_pReal, 0.0_pReal, 0.0_pReal] ! original axis is now +x-axis
case('-x')
texture_transformation(j,1:3,t) = [-1.0_pReal, 0.0_pReal, 0.0_pReal] ! original axis is now -x-axis
case('y', '+y')
texture_transformation(j,1:3,t) = [ 0.0_pReal, 1.0_pReal, 0.0_pReal] ! original axis is now +y-axis
case('-y')
texture_transformation(j,1:3,t) = [ 0.0_pReal,-1.0_pReal, 0.0_pReal] ! original axis is now -y-axis
case('z', '+z')
texture_transformation(j,1:3,t) = [ 0.0_pReal, 0.0_pReal, 1.0_pReal] ! original axis is now +z-axis
case('-z')
texture_transformation(j,1:3,t) = [ 0.0_pReal, 0.0_pReal,-1.0_pReal] ! original axis is now -z-axis
case default
call IO_error(157,t)
end select
enddo
if(dNeq(math_det33(texture_transformation(1:3,1:3,t)),1.0_pReal)) call IO_error(157,t)
endif
if (config_texture(t)%keyExists('(gauss)')) then
gauss = gauss + 1
strings = config_texture(t)%getStrings('(gauss)',raw= .true.)
do i = 1 , size(strings)
chunkPos = IO_stringPos(strings(i))
do j = 1,9,2
select case (IO_stringValue(strings(i),chunkPos,j))
case('phi1')
texture_Gauss(1,gauss,t) = IO_floatValue(strings(i),chunkPos,j+1)*inRad
case('phi')
texture_Gauss(2,gauss,t) = IO_floatValue(strings(i),chunkPos,j+1)*inRad
case('phi2')
texture_Gauss(3,gauss,t) = IO_floatValue(strings(i),chunkPos,j+1)*inRad
end select
do t=1, size(config_texture)
strings = config_texture(t)%getStrings('(gauss)',raw= .true.)
do i = 1 , size(strings)
chunkPos = IO_stringPos(strings(i))
do j = 1,9,2
select case (IO_stringValue(strings(i),chunkPos,j))
case('phi1')
texture_Gauss(1,t) = IO_floatValue(strings(i),chunkPos,j+1)*inRad
case('phi')
texture_Gauss(2,t) = IO_floatValue(strings(i),chunkPos,j+1)*inRad
case('phi2')
texture_Gauss(3,t) = IO_floatValue(strings(i),chunkPos,j+1)*inRad
end select
enddo
enddo
endif
enddo
call config_deallocate('material.config/texture')
enddo
if (config_texture(t)%keyExists('axes')) then
strings = config_texture(t)%getStrings('axes')
do j = 1, 3 ! look for "x", "y", and "z" entries
select case (strings(j))
case('x', '+x')
texture_transformation(j,1:3) = [ 1.0_pReal, 0.0_pReal, 0.0_pReal] ! original axis is now +x-axis
case('-x')
texture_transformation(j,1:3) = [-1.0_pReal, 0.0_pReal, 0.0_pReal] ! original axis is now -x-axis
case('y', '+y')
texture_transformation(j,1:3) = [ 0.0_pReal, 1.0_pReal, 0.0_pReal] ! original axis is now +y-axis
case('-y')
texture_transformation(j,1:3) = [ 0.0_pReal,-1.0_pReal, 0.0_pReal] ! original axis is now -y-axis
case('z', '+z')
texture_transformation(j,1:3) = [ 0.0_pReal, 0.0_pReal, 1.0_pReal] ! original axis is now +z-axis
case('-z')
texture_transformation(j,1:3) = [ 0.0_pReal, 0.0_pReal,-1.0_pReal] ! original axis is now -z-axis
case default
call IO_error(157,t)
end select
enddo
if(dNeq(math_det33(texture_transformation),1.0_pReal)) call IO_error(157,t)
call eulers%fromEulerAngles(texture_Gauss(:,t))
texture_Gauss(:,t) = math_RtoEuler(matmul(eulers%asRotationMatrix(),texture_transformation))
endif
enddo
call config_deallocate('material.config/texture')
end subroutine material_parseTexture
@ -876,25 +873,16 @@ subroutine material_populateGrains
homog = discretization_homogenizationAt(e)
micro = discretization_microstructureAt(e)
do c = 1, homogenization_Ngrains(homog)
material_phase(c,i,e) = microstructure_phase(c,micro)
material_texture(c,i,e) = microstructure_texture(c,micro)
material_EulerAngles(1:3,c,i,e) = texture_Gauss(1:3,1,material_texture(c,i,e))
material_EulerAngles(1:3,c,i,e) = math_RtoEuler( & ! translate back to Euler angles
matmul( & ! pre-multiply
math_EulertoR(material_EulerAngles(1:3,c,i,e)), & ! face-value orientation
texture_transformation(1:3,1:3,material_texture(c,i,e)) & ! and transformation matrix
) &
)
material_phase(c,i,e) = microstructure_phase(c,micro)
material_texture(c,i,e) = microstructure_texture(c,micro)
material_EulerAngles(1:3,c,i,e) = texture_Gauss(1:3,material_texture(c,i,e))
enddo
enddo
enddo
deallocate(texture_transformation)
deallocate(microstructure_phase)
deallocate(microstructure_texture)
call config_deallocate('material.config/microstructure')
end subroutine material_populateGrains
end module material

View File

@ -872,25 +872,27 @@ subroutine calcCells(thisMesh,elem,connectivity_elem)
class(tMesh) :: thisMesh
type(tElement) :: elem
integer(pInt),dimension(:,:), intent(in) :: connectivity_elem
integer(pInt),dimension(:,:), allocatable :: con_elem,temp,con,parentsAndWeights,candidates_global
integer(pInt),dimension(:), allocatable :: l, nodes, candidates_local
integer(pInt),dimension(:,:,:), allocatable :: con_cell,connectivity_cell
integer(pInt),dimension(:,:), allocatable :: sorted,test,connectivity_cell_reshape
integer,dimension(:,:), intent(in) :: connectivity_elem
integer,dimension(:,:), allocatable :: con_elem,temp,con,parentsAndWeights,candidates_global
integer,dimension(:), allocatable :: l, nodes, candidates_local
integer,dimension(:,:,:), allocatable :: con_cell,connectivity_cell
integer,dimension(:,:), allocatable :: sorted,test,connectivity_cell_reshape
real(pReal), dimension(:,:), allocatable :: coordinates,nodes5
integer(pInt) :: e, n, c, p, s,u,i,m,j,nParentNodes,nCellNode,ierr
integer :: e, n, c, p, s,u,i,m,j,nParentNodes,nCellNode,ierr,Nelem,candidateID
Nelem = thisMesh%Nelems
!---------------------------------------------------------------------------------------------------
! initialize global connectivity to negative local connectivity
allocate(connectivity_cell(thisMesh%elem%NcellNodesPerCell,thisMesh%elem%nIPs,thisMesh%Nelems))
connectivity_cell = -spread(thisMesh%elem%cell,3,thisMesh%Nelems) ! local cell node ID
allocate(connectivity_cell(elem%NcellNodesPerCell,elem%nIPs,Nelem))
connectivity_cell = -spread(elem%cell,3,Nelem) ! local cell node ID
!---------------------------------------------------------------------------------------------------
! set connectivity of cell nodes that conincide with FE nodes (defined by 1 parent node)
! change to global node ID
do e = 1, thisMesh%Nelems
do c = 1, thisMesh%elem%NcellNodes
realNode: if (count(thisMesh%elem%cellNodeParentNodeWeights(:,c) /= 0_pInt) == 1_pInt) then
! set connectivity of cell nodes that coincide with FE nodes (defined by 1 parent node)
! and renumber local (negative) to global (positive) node ID
do e = 1, Nelem
do c = 1, elem%NcellNodes
realNode: if (count(elem%cellNodeParentNodeWeights(:,c) /= 0) == 1) then
where(connectivity_cell(:,:,e) == -c)
connectivity_cell(:,:,e) = connectivity_elem(c,e)
end where
@ -900,109 +902,114 @@ subroutine calcCells(thisMesh,elem,connectivity_elem)
nCellNode = thisMesh%nNodes
do nParentNodes = 2, thisMesh%elem%nNodes
!---------------------------------------------------------------------------------------------------
! set connectivity of cell nodes that are defined by 2,...,nNodes real nodes
do nParentNodes = 2, elem%nNodes
! get IDs of local cell nodes that are defined by the current number of parent nodes
candidates_local = [integer(pInt)::]
do c = 1, thisMesh%elem%NcellNodes
if (count(thisMesh%elem%cellNodeParentNodeWeights(:,c) /= 0_pInt) == nParentNodes) &
candidates_local = [candidates_local,c]
enddo
s = size(candidates_local)
if (allocated(candidates_global)) deallocate(candidates_global)
allocate(candidates_global(nParentNodes*2_pInt+2_pInt,s*thisMesh%Nelems))
parentsAndWeights = reshape([(0_pInt, i = 1_pInt,2_pInt*nParentNodes)],[nParentNodes,2])
do e = 1_pInt, thisMesh%Nelems
do i = 1_pInt, s
c = candidates_local(i)
m = 0_pInt
do p = 1_pInt, size(thisMesh%elem%cellNodeParentNodeWeights(:,c))
if (thisMesh%elem%cellNodeParentNodeWeights(p,c) /= 0_pInt) then ! real node 'c' partly defines cell node 'n'
m = m + 1_pInt
parentsAndWeights(m,1:2) = [connectivity_elem(p,e),thisMesh%elem%cellNodeParentNodeWeights(p,c)]
endif
enddo
! store (and order) real nodes and their weights together with the element number and local ID
do p = 1_pInt, nParentNodes
m = maxloc(parentsAndWeights(:,1),1)
candidates_global(p, (e-1)*s+i) = parentsAndWeights(m,1)
parentsAndWeights(m,1) = -huge(i) ! out of the competition
candidates_global(p+nParentNodes,(e-1)*s+i) = parentsAndWeights(m,2)
candidates_global(nParentNodes*2+1:nParentNodes*2+2,(e-1)*s+i) = [e,c]
enddo
enddo
enddo
! sort according to real node IDs (from left to right)
call math_sort(candidates_global,sortDim=1)
do p = 2, nParentNodes-1 ! why -1?
n = 1
do while(n <= s*thisMesh%Nelems)
j=0
do while (n+j<= s*thisMesh%Nelems)
if (candidates_global(p-1,n+j)/=candidates_global(p-1,n)) exit
j = j + 1
! get IDs of local cell nodes that are defined by the current number of parent nodes
candidates_local = [integer::]
do c = 1, elem%NcellNodes
if (count(elem%cellNodeParentNodeWeights(:,c) /= 0) == nParentNodes) &
candidates_local = [candidates_local,c]
enddo
s = size(candidates_local)
if (allocated(candidates_global)) deallocate(candidates_global)
allocate(candidates_global(nParentNodes*2+2,s*Nelem)) ! stores parent node ID + weight together with element ID and cellnode id (local)
parentsAndWeights = reshape([(0, i = 1,2*nParentNodes)],[nParentNodes,2]) ! allocate without deallocate
do e = 1, Nelem
do i = 1, size(candidates_local)
candidateID = (e-1)*size(candidates_local)+i ! including duplicates (Nelem*size(candidates_local))
c = candidates_local(i) ! c is local cellnode ID for connectivity
p = 0
do j = 1, size(elem%cellNodeParentNodeWeights(:,c))
if (elem%cellNodeParentNodeWeights(j,c) /= 0) then ! real node 'j' partly defines cell node 'c'
p = p + 1
parentsAndWeights(p,1:2) = [connectivity_elem(j,e),elem%cellNodeParentNodeWeights(j,c)]
endif
enddo
! store (and order) real node IDs and their weights together with the element number and local ID
do p = 1, nParentNodes
m = maxloc(parentsAndWeights(:,1),1)
candidates_global(p, candidateID) = parentsAndWeights(m,1)
candidates_global(p+nParentNodes, candidateID) = parentsAndWeights(m,2)
candidates_global(nParentNodes*2+1:nParentNodes*2+2,candidateID) = [e,c]
parentsAndWeights(m,1) = -huge(parentsAndWeights(m,1)) ! out of the competition
enddo
enddo
e = n+j-1
if (any(candidates_global(p,n:e)/=candidates_global(p,n))) &
call math_sort(candidates_global(:,n:e),sortDim=p)
n = e+1
enddo
enddo
! find duplicates (trivial for sorted IDs)
i = 0
n = 1
do while(n <= s*thisMesh%Nelems)
enddo
! sort according to real node IDs + weight (from left to right)
call math_sort(candidates_global,sortDim=1) ! sort according to first column
do p = 2, nParentNodes*2
n = 1
do while(n <= size(candidates_local)*Nelem)
j=0
do while (n+j<= size(candidates_local)*Nelem)
if (candidates_global(p-1,n+j)/=candidates_global(p-1,n)) exit
j = j + 1
enddo
e = n+j-1
if (any(candidates_global(p,n:e)/=candidates_global(p,n))) &
call math_sort(candidates_global(:,n:e),sortDim=p)
n = e+1
enddo
enddo
! count unique cell nodes (trivial for sorted IDs + weights)
i = 0 ! counts unique cell nodes (defined by current number of parents)
n = 1
do while(n <= size(candidates_local)*Nelem)
j=0
do while (n+j<= s*thisMesh%Nelems)
do while (n+j<= size(candidates_local)*Nelem)
if (any(candidates_global(1:2*nParentNodes,n+j)/=candidates_global(1:2*nParentNodes,n))) exit
j = j + 1
enddo
i=i+1
i = i+1
n = n+j
enddo
! calculate coordinates of cell nodes and insert their ID into the cell conectivity
coordinates = reshape([(0.0_pReal,j = 1, 3*i)], [3,i])
i = 1
n = 1
do while(n <= size(candidates_local)*Nelem)
j=0
parentsAndWeights(:,1) = candidates_global(1:nParentNodes,n+j)
parentsAndWeights(:,2) = candidates_global(nParentNodes+1:nParentNodes*2,n+j)
e = candidates_global(nParentNodes*2+1,n+j)
c = candidates_global(nParentNodes*2+2,n+j)
do m = 1, nParentNodes
coordinates(:,i) = coordinates(:,i) &
+ thisMesh%node_0(:,parentsAndWeights(m,1)) * real(parentsAndWeights(m,2),pReal)
enddo
p = i ! ToDo: Hack
! calculate coordinates of cell nodes and insert their ID into the cell conectivity
coordinates = reshape([(0.0_pReal,i = 1, 3*s*thisMesh%Nelems)], [3,i])
i = 0
n = 1
do while(n <= s*thisMesh%Nelems)
j=0
parentsAndWeights(:,1) = candidates_global(1:nParentNodes,n+j)
parentsAndWeights(:,2) = candidates_global(nParentNodes+1:nParentNodes*2,n+j)
e = candidates_global(nParentNodes*2+1,n+j)
c = candidates_global(nParentNodes*2+2,n+j)
do m = 1, nParentNodes
coordinates(:,i+1) = coordinates(:,i+1) &
+ thisMesh%node_0(:,parentsAndWeights(m,1)) * real(parentsAndWeights(m,2),pReal)
enddo
coordinates(:,i+1) = coordinates(:,i+1)/real(sum(parentsAndWeights(:,2)),pReal)
do while (n+j<= s*thisMesh%Nelems)
if (any(candidates_global(1:2*nParentNodes,n+j)/=candidates_global(1:2*nParentNodes,n))) exit
where (connectivity_cell(:,:,candidates_global(nParentNodes*2+1,n+j)) == -candidates_global(nParentNodes*2+2,n+j))
connectivity_cell(:,:,candidates_global(nParentNodes*2+1,n+j)) = i+1+nCellNode
end where
j = j + 1
enddo
i=i+1
n = n+j
coordinates(:,i) = coordinates(:,i)/real(sum(parentsAndWeights(:,2)),pReal)
do while (n+j<= size(candidates_local)*Nelem)
if (any(candidates_global(1:2*nParentNodes,n+j)/=candidates_global(1:2*nParentNodes,n))) exit
where (connectivity_cell(:,:,candidates_global(nParentNodes*2+1,n+j)) == -candidates_global(nParentNodes*2+2,n+j)) ! still locally defined
connectivity_cell(:,:,candidates_global(nParentNodes*2+1,n+j)) = nCellNode + i ! gets current new cell node id
end where
j = j + 1
enddo
nCellNode = nCellNode + p
if (p/=0) nodes5 = reshape([nodes5,coordinates],[3,nCellNode])
enddo
thisMesh%node_0 = nodes5
mesh_cell2 = connectivity_cell
connectivity_cell_reshape = reshape(connectivity_cell,[thisMesh%elem%NcellNodesPerCell,thisMesh%elem%nIPs*thisMesh%Nelems])
i=i+1
n = n+j
enddo
nCellNode = nCellNode + i
if (i/=0) nodes5 = reshape([nodes5,coordinates],[3,nCellNode])
enddo
thisMesh%node_0 = nodes5
mesh_cell2 = connectivity_cell
connectivity_cell_reshape = reshape(connectivity_cell,[elem%NcellNodesPerCell,elem%nIPs*thisMesh%Nelems])
#if defined(DAMASK_HDF5)
call results_openJobFile

View File

@ -65,6 +65,7 @@ module rotations
procedure, public :: asRotationMatrix
!------------------------------------------
procedure, public :: fromRotationMatrix
procedure, public :: fromEulerAngles
!------------------------------------------
procedure, public :: rotVector
procedure, public :: rotTensor
@ -143,7 +144,16 @@ subroutine fromRotationMatrix(self,om)
self%q = om2qu(om)
end subroutine
!---------------------------------------------------------------------------------------------------
subroutine fromEulerAngles(self,eu)
class(rotation), intent(out) :: self
real(pReal), dimension(3), intent(in) :: eu
self%q = eu2qu(eu)
end subroutine
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University