223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

initialize mesh and element

This commit is contained in:
Martin Diehl 2019-02-01 12:24:23 +01:00
parent b87a09a466
commit 5f8b110f63
4 changed files with 70 additions and 162 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

109
src/mesh_abaqus.f90
View File

@ -389,7 +389,7 @@ integer(pInt), dimension(:,:), allocatable, private :: &
mesh_abaqus_build_elements mesh_abaqus_build_elements
type, public, extends(tMesh) :: tMesh_Abaqus type, public, extends(tMesh) :: tMesh_abaqus
integer(pInt):: & integer(pInt):: &
mesh_Nelems, & !< total number of elements in mesh (including non-DAMASK elements) mesh_Nelems, & !< total number of elements in mesh (including non-DAMASK elements)
@ -406,16 +406,22 @@ integer(pInt), dimension(:,:), allocatable, private :: &
logical:: noPart !< for cases where the ABAQUS input file does not use part/assembly information logical:: noPart !< for cases where the ABAQUS input file does not use part/assembly information
contains contains
procedure :: init=>tMesh_abaqus_init procedure, pass(self) :: tMesh_abaqus_init
end type tMesh_Abaqus generic, public :: init => tMesh_abaqus_init
end type tMesh_abaqus
type(tMesh_Abaqus), public, protected :: theMesh type(tMesh_abaqus), public, protected :: theMesh
contains contains
subroutine tMesh_abaqus_init(self) subroutine tMesh_abaqus_init(self,elemType,nodes)
implicit none implicit none
class(tMesh_abaqus) :: self class(tMesh_abaqus) :: self
real(pReal), dimension(:,:), intent(in) :: nodes
integer(pInt), intent(in) :: elemType
call self%tMesh%init('mesh',elemType,nodes)
end subroutine tMesh_abaqus_init end subroutine tMesh_abaqus_init
@ -537,7 +543,7 @@ subroutine mesh_init(ip,el)
mesh_microstructureAt = mesh_element(4,:) mesh_microstructureAt = mesh_element(4,:)
mesh_CPnodeID = mesh_element(5:4+mesh_NipsPerElem,:) mesh_CPnodeID = mesh_element(5:4+mesh_NipsPerElem,:)
!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!
call theMesh%init(mesh_element(2,1),mesh_node0)
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief check if the input file for Abaqus contains part info !> @brief check if the input file for Abaqus contains part info
@ -859,97 +865,6 @@ pure function mesh_cellCenterCoordinates(ip,el)
end function mesh_cellCenterCoordinates end function mesh_cellCenterCoordinates
!--------------------------------------------------------------------------------------------------
!> @brief builds mesh of (distorted) cubes for given coordinates (= center of the cubes)
!--------------------------------------------------------------------------------------------------
function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes)
use debug, only: &
debug_mesh, &
debug_level, &
debug_levelBasic
use math, only: &
math_mul33x3
implicit none
real(pReal), intent(in), dimension(:,:,:,:) :: &
centres
real(pReal), dimension(3,size(centres,2)+1,size(centres,3)+1,size(centres,4)+1) :: &
nodes
real(pReal), intent(in), dimension(3) :: &
gDim
real(pReal), intent(in), dimension(3,3) :: &
Favg
real(pReal), dimension(3,size(centres,2)+2,size(centres,3)+2,size(centres,4)+2) :: &
wrappedCentres
integer(pInt) :: &
i,j,k,n
integer(pInt), dimension(3), parameter :: &
diag = 1_pInt
integer(pInt), dimension(3) :: &
shift = 0_pInt, &
lookup = 0_pInt, &
me = 0_pInt, &
iRes = 0_pInt
integer(pInt), dimension(3,8) :: &
neighbor = reshape([ &
0_pInt, 0_pInt, 0_pInt, &
1_pInt, 0_pInt, 0_pInt, &
1_pInt, 1_pInt, 0_pInt, &
0_pInt, 1_pInt, 0_pInt, &
0_pInt, 0_pInt, 1_pInt, &
1_pInt, 0_pInt, 1_pInt, &
1_pInt, 1_pInt, 1_pInt, &
0_pInt, 1_pInt, 1_pInt ], [3,8])
!--------------------------------------------------------------------------------------------------
! initializing variables
iRes = [size(centres,2),size(centres,3),size(centres,4)]
nodes = 0.0_pReal
wrappedCentres = 0.0_pReal
!--------------------------------------------------------------------------------------------------
! report
if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then
write(6,'(a)') ' Meshing cubes around centroids'
write(6,'(a,3(e12.5))') ' Dimension: ', gDim
write(6,'(a,3(i5))') ' Resolution:', iRes
endif
!--------------------------------------------------------------------------------------------------
! building wrappedCentres = centroids + ghosts
wrappedCentres(1:3,2_pInt:iRes(1)+1_pInt,2_pInt:iRes(2)+1_pInt,2_pInt:iRes(3)+1_pInt) = centres
do k = 0_pInt,iRes(3)+1_pInt
do j = 0_pInt,iRes(2)+1_pInt
do i = 0_pInt,iRes(1)+1_pInt
if (k==0_pInt .or. k==iRes(3)+1_pInt .or. & ! z skin
j==0_pInt .or. j==iRes(2)+1_pInt .or. & ! y skin
i==0_pInt .or. i==iRes(1)+1_pInt ) then ! x skin
me = [i,j,k] ! me on skin
shift = sign(abs(iRes+diag-2_pInt*me)/(iRes+diag),iRes+diag-2_pInt*me)
lookup = me-diag+shift*iRes
wrappedCentres(1:3,i+1_pInt, j+1_pInt, k+1_pInt) = &
centres(1:3,lookup(1)+1_pInt,lookup(2)+1_pInt,lookup(3)+1_pInt) &
- math_mul33x3(Favg, real(shift,pReal)*gDim)
endif
enddo; enddo; enddo
!--------------------------------------------------------------------------------------------------
! averaging
do k = 0_pInt,iRes(3); do j = 0_pInt,iRes(2); do i = 0_pInt,iRes(1)
do n = 1_pInt,8_pInt
nodes(1:3,i+1_pInt,j+1_pInt,k+1_pInt) = &
nodes(1:3,i+1_pInt,j+1_pInt,k+1_pInt) + wrappedCentres(1:3,i+1_pInt+neighbor(1,n), &
j+1_pInt+neighbor(2,n), &
k+1_pInt+neighbor(3,n) )
enddo
enddo; enddo; enddo
nodes = nodes/8.0_pReal
end function mesh_nodesAroundCentres
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Count overall number of nodes and elements in mesh and stores them in !> @brief Count overall number of nodes and elements in mesh and stores them in
!! 'mesh_Nelems' and 'mesh_Nnodes' !! 'mesh_Nelems' and 'mesh_Nnodes'

20
src/mesh_base.f90
View File

@ -29,7 +29,7 @@ module mesh_base
node !< node x,y,z coordinates (deformed) node !< node x,y,z coordinates (deformed)
integer(pInt), dimension(:,:), allocatable, public :: & integer(pInt), dimension(:,:), allocatable, public :: &
cellnodeParent !< cellnode's parent element ID, cellnode's intra-element ID cellnodeParent !< cellnode's parent element ID, cellnode's intra-element ID
character(pStringLen) :: solver = "undefined" character(pStringLen) :: type = "n/a"
integer(pInt) :: & integer(pInt) :: &
Nnodes, & !< total number of nodes in mesh Nnodes, & !< total number of nodes in mesh
Nelems = -1_pInt, & Nelems = -1_pInt, &
@ -43,6 +43,24 @@ module mesh_base
microstructureAt microstructureAt
integer(pInt), dimension(:,:), allocatable, public :: & integer(pInt), dimension(:,:), allocatable, public :: &
connectivity connectivity
contains
procedure, pass(self) :: tMesh_base_init
generic, public :: init => tMesh_base_init
end type tMesh end type tMesh
contains
subroutine tMesh_base_init(self,meshType,elemType,nodes)
implicit none
class(tMesh) :: self
character(len=*), intent(in) :: meshType
integer(pInt), intent(in) :: elemType
real(pReal), dimension(:,:), intent(in) :: nodes
self%type = meshType
call self%elem%init(elemType)
self%node0 = nodes
end subroutine tMesh_base_init
end module mesh_base end module mesh_base

86
src/mesh_grid.f90
View File

@ -270,14 +270,11 @@ integer(pInt), dimension(:,:), allocatable, private :: &
mesh_init, & mesh_init, &
mesh_cellCenterCoordinates mesh_cellCenterCoordinates
private :: & private :: &
mesh_build_cellconnectivity, & mesh_build_cellconnectivity, &
mesh_build_ipAreas, & mesh_build_ipAreas, &
mesh_build_FEdata, & mesh_build_FEdata, &
mesh_spectral_getHomogenization, & mesh_spectral_getHomogenization, &
mesh_spectral_count, &
mesh_spectral_count_cpSizes, &
mesh_spectral_build_nodes, & mesh_spectral_build_nodes, &
mesh_spectral_build_elements, & mesh_spectral_build_elements, &
mesh_spectral_build_ipNeighborhood, & mesh_spectral_build_ipNeighborhood, &
@ -302,19 +299,21 @@ integer(pInt), dimension(:,:), allocatable, private :: &
size3offset size3offset
contains contains
procedure :: init => tMesh_grid_init procedure, pass(self) :: tMesh_grid_init
generic, public :: init => tMesh_grid_init
end type tMesh_grid end type tMesh_grid
type(tMesh_grid), public, protected :: theMesh type(tMesh_grid), public, protected :: theMesh
contains contains
subroutine tMesh_grid_init(self) subroutine tMesh_grid_init(self,nodes)
implicit none implicit none
class(tMesh_grid) :: self class(tMesh_grid) :: self
real(pReal), dimension(:,:), intent(in) :: nodes
call self%elem%init(10_pInt) call self%tMesh%init('grid',10_pInt,nodes)
end subroutine tMesh_grid_init end subroutine tMesh_grid_init
@ -364,7 +363,8 @@ subroutine mesh_init(ip,el)
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
call theMesh%init
call mesh_build_FEdata ! get properties of the different types of elements call mesh_build_FEdata ! get properties of the different types of elements
mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh
@ -389,13 +389,23 @@ subroutine mesh_init(ip,el)
grid3Offset = int(local_K_offset,pInt) grid3Offset = int(local_K_offset,pInt)
size3 = geomSize(3)*real(grid3,pReal) /real(grid(3),pReal) size3 = geomSize(3)*real(grid3,pReal) /real(grid(3),pReal)
size3Offset = geomSize(3)*real(grid3Offset,pReal)/real(grid(3),pReal) size3Offset = geomSize(3)*real(grid3Offset,pReal)/real(grid(3),pReal)
if (myDebug) write(6,'(a)') ' Grid partitioned'; flush(6) mesh_NcpElems= product(grid(1:2))*grid3
call mesh_spectral_count() mesh_NcpElemsGlobal = product(grid)
if (myDebug) write(6,'(a)') ' Counted nodes/elements'; flush(6)
call mesh_spectral_count_cpSizes mesh_Nnodes = product(grid(1:2) + 1_pInt)*(grid3 + 1_pInt)
if (myDebug) write(6,'(a)') ' Built CP statistics'; flush(6)
call mesh_spectral_build_nodes() call mesh_spectral_build_nodes()
if (myDebug) write(6,'(a)') ' Built nodes'; flush(6) if (myDebug) write(6,'(a)') ' Built nodes'; flush(6)
call theMesh%init(mesh_node)
! For compatibility
mesh_maxNips = theMesh%elem%nIPs
mesh_maxNipNeighbors = theMesh%elem%nIPneighbors
mesh_maxNcellnodes = theMesh%elem%Ncellnodes
call mesh_spectral_build_elements(FILEUNIT) call mesh_spectral_build_elements(FILEUNIT)
if (myDebug) write(6,'(a)') ' Built elements'; flush(6) if (myDebug) write(6,'(a)') ' Built elements'; flush(6)
call mesh_build_cellconnectivity call mesh_build_cellconnectivity
@ -434,8 +444,6 @@ subroutine mesh_init(ip,el)
mesh_microstructureAt = mesh_element(4,:) mesh_microstructureAt = mesh_element(4,:)
mesh_CPnodeID = mesh_element(5:4+mesh_NipsPerElem,:) mesh_CPnodeID = mesh_element(5:4+mesh_NipsPerElem,:)
!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!
end subroutine mesh_init end subroutine mesh_init
@ -563,10 +571,9 @@ subroutine mesh_build_ipVolumes
integer(pInt) :: e,t,g,c,i,m,f,n integer(pInt) :: e,t,g,c,i,m,f,n
real(pReal), dimension(FE_maxNcellnodesPerCellface,FE_maxNcellfaces) :: subvolume real(pReal), dimension(FE_maxNcellnodesPerCellface,FE_maxNcellfaces) :: subvolume
if (.not. allocated(mesh_ipVolume)) then
allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems)) allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems),source=0.0_pReal)
mesh_ipVolume = 0.0_pReal
endif
!$OMP PARALLEL DO PRIVATE(t,g,c,m,subvolume) !$OMP PARALLEL DO PRIVATE(t,g,c,m,subvolume)
do e = 1_pInt,mesh_NcpElems ! loop over cpElems do e = 1_pInt,mesh_NcpElems ! loop over cpElems
@ -894,43 +901,6 @@ integer(pInt) function mesh_spectral_getHomogenization(fileUnit)
end function mesh_spectral_getHomogenization end function mesh_spectral_getHomogenization
!--------------------------------------------------------------------------------------------------
!> @brief Count overall number of nodes and elements in mesh and stores them in
!! 'mesh_Nelems', 'mesh_Nnodes' and 'mesh_NcpElems'
!--------------------------------------------------------------------------------------------------
subroutine mesh_spectral_count()
implicit none
mesh_NcpElems= product(grid(1:2))*grid3
mesh_Nnodes = product(grid(1:2) + 1_pInt)*(grid3 + 1_pInt)
mesh_NcpElemsGlobal = product(grid)
end subroutine mesh_spectral_count
!--------------------------------------------------------------------------------------------------
!> @brief Gets maximum count of nodes, IPs, IP neighbors, and subNodes among cpElements.
!! Sets global values 'mesh_maxNips', 'mesh_maxNipNeighbors',
!! and 'mesh_maxNcellnodes'
!--------------------------------------------------------------------------------------------------
subroutine mesh_spectral_count_cpSizes
implicit none
integer(pInt) :: t,g,c
t = 10_pInt
g = FE_geomtype(t)
c = FE_celltype(g)
mesh_maxNips = FE_Nips(g)
mesh_maxNipNeighbors = FE_NipNeighbors(c)
mesh_maxNcellnodes = FE_Ncellnodes(g)
end subroutine mesh_spectral_count_cpSizes
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Store x,y,z coordinates of all nodes in mesh. !> @brief Store x,y,z coordinates of all nodes in mesh.
!! Allocates global arrays 'mesh_node0' and 'mesh_node' !! Allocates global arrays 'mesh_node0' and 'mesh_node'
@ -941,7 +911,6 @@ subroutine mesh_spectral_build_nodes()
integer(pInt) :: n integer(pInt) :: n
allocate (mesh_node0 (3,mesh_Nnodes), source = 0.0_pReal) allocate (mesh_node0 (3,mesh_Nnodes), source = 0.0_pReal)
allocate (mesh_node (3,mesh_Nnodes), source = 0.0_pReal)
forall (n = 0_pInt:mesh_Nnodes-1_pInt) forall (n = 0_pInt:mesh_Nnodes-1_pInt)
mesh_node0(1,n+1_pInt) = mesh_unitlength * & mesh_node0(1,n+1_pInt) = mesh_unitlength * &
@ -986,7 +955,6 @@ subroutine mesh_spectral_build_elements(fileUnit)
headerLength = 0_pInt, & headerLength = 0_pInt, &
maxDataPerLine, & maxDataPerLine, &
homog, & homog, &
elemType, &
elemOffset elemOffset
integer(pInt), dimension(:), allocatable :: & integer(pInt), dimension(:), allocatable :: &
microstructures, & microstructures, &
@ -1047,13 +1015,13 @@ subroutine mesh_spectral_build_elements(fileUnit)
enddo enddo
enddo enddo
elemType = 10_pInt
elemOffset = product(grid(1:2))*grid3Offset elemOffset = product(grid(1:2))*grid3Offset
e = 0_pInt e = 0_pInt
do while (e < mesh_NcpElems) ! fill expected number of elements, stop at end of data (or blank line!) do while (e < mesh_NcpElems) ! fill expected number of elements, stop at end of data (or blank line!)
e = e+1_pInt ! valid element entry e = e+1_pInt ! valid element entry
mesh_element( 1,e) = -1_pInt ! DEPRECATED mesh_element( 1,e) = -1_pInt ! DEPRECATED
mesh_element( 2,e) = elemType ! elem type mesh_element( 2,e) = 10_pInt
mesh_element( 3,e) = homog ! homogenization mesh_element( 3,e) = homog ! homogenization
mesh_element( 4,e) = microGlobal(e+elemOffset) ! microstructure mesh_element( 4,e) = microGlobal(e+elemOffset) ! microstructure
mesh_element( 5,e) = e + (e-1_pInt)/grid(1) + & mesh_element( 5,e) = e + (e-1_pInt)/grid(1) + &

17
src/mesh_marc.f90
View File

@ -417,8 +417,9 @@ type, public, extends(tMesh) :: tMesh_marc
integer(pInt), dimension(:,:), allocatable :: & integer(pInt), dimension(:,:), allocatable :: &
mesh_mapElemSet !< list of elements in elementSet mesh_mapElemSet !< list of elements in elementSet
contains contains
procedure :: init => tMesh_marc_init procedure, pass(self) :: tMesh_marc_init
generic, public :: init => tMesh_marc_init
end type tMesh_marc end type tMesh_marc
type(tMesh_marc), public, protected :: theMesh type(tMesh_marc), public, protected :: theMesh
@ -426,10 +427,15 @@ end type tMesh_marc
contains contains
subroutine tMesh_marc_init(self) subroutine tMesh_marc_init(self,elemType,nodes)
implicit none implicit none
class(tMesh_marc) :: self class(tMesh_marc) :: self
real(pReal), dimension(:,:), intent(in) :: nodes
integer(pInt), intent(in) :: elemType
call self%tMesh%init('mesh',elemType,nodes)
end subroutine tMesh_marc_init end subroutine tMesh_marc_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -553,7 +559,8 @@ subroutine mesh_init(ip,el)
mesh_microstructureAt = mesh_element(4,:) mesh_microstructureAt = mesh_element(4,:)
mesh_CPnodeID = mesh_element(5:4+mesh_NipsPerElem,:) mesh_CPnodeID = mesh_element(5:4+mesh_NipsPerElem,:)
!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!
call theMesh%init(mesh_element(2,1),mesh_node0)
end subroutine mesh_init end subroutine mesh_init