diff --git a/src/material.f90 b/src/material.f90 index e63f05cd3..32a0f335d 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -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 diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 254fee5bc..97ba2634d 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -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 diff --git a/src/rotations.f90 b/src/rotations.f90 index 3a64f27b9..5090261e6 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -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