cleaning
This commit is contained in:
Martin Diehl
parent
e3b16639bf
commit
e0cb1a87cd
|
|
@ -125,10 +125,6 @@ subroutine mesh_init(ip,el)
|
||||||
integer, dimension(:), allocatable :: &
|
integer, dimension(:), allocatable :: &
|
||||||
marc_matNumber !< array of material numbers for hypoelastic material (Marc only)
|
marc_matNumber !< array of material numbers for hypoelastic material (Marc only)
|
||||||
logical :: myDebug
|
logical :: myDebug
|
||||||
real(pReal), dimension(:,:,:), allocatable:: &
|
|
||||||
mesh_ipArea !< area of interface to neighboring IP (initially!)
|
|
||||||
real(pReal),dimension(:,:,:,:), allocatable :: &
|
|
||||||
mesh_ipAreaNormal !< area normal of interface to neighboring IP (initially!)
|
|
||||||
|
|
||||||
real(pReal), dimension(:,:), allocatable :: &
|
real(pReal), dimension(:,:), allocatable :: &
|
||||||
ip_reshaped
|
ip_reshaped
|
||||||
|
|
@ -195,9 +191,6 @@ subroutine mesh_init(ip,el)
|
||||||
call mesh_build_ipCoordinates
|
call mesh_build_ipCoordinates
|
||||||
if (myDebug) write(6,'(a)') ' Built IP coordinates'; flush(6)
|
if (myDebug) write(6,'(a)') ' Built IP coordinates'; flush(6)
|
||||||
|
|
||||||
allocate(mesh_ipArea(theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems))
|
|
||||||
allocate(mesh_ipAreaNormal(3,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems))
|
|
||||||
call mesh_build_ipAreas(mesh_ipArea,mesh_ipAreaNormal)
|
|
||||||
if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6)
|
if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6)
|
||||||
|
|
||||||
call IP_neighborhood2
|
call IP_neighborhood2
|
||||||
|
|
@ -230,10 +223,7 @@ subroutine mesh_init(ip,el)
|
||||||
'nodal coordinates','m')
|
'nodal coordinates','m')
|
||||||
call results_closeJobFile()
|
call results_closeJobFile()
|
||||||
#endif
|
#endif
|
||||||
call geometry_plastic_nonlocal_setIPvolume(IPvolume())
|
|
||||||
call geometry_plastic_nonlocal_setIPneighborhood(mesh_ipNeighborhood2)
|
call geometry_plastic_nonlocal_setIPneighborhood(mesh_ipNeighborhood2)
|
||||||
call geometry_plastic_nonlocal_setIParea(mesh_ipArea)
|
|
||||||
call geometry_plastic_nonlocal_setIPareaNormal(mesh_ipAreaNormal)
|
|
||||||
|
|
||||||
end subroutine mesh_init
|
end subroutine mesh_init
|
||||||
|
|
||||||
|
|
@ -1018,68 +1008,6 @@ function mesh_build_cellnodes()
|
||||||
end function mesh_build_cellnodes
|
end function mesh_build_cellnodes
|
||||||
|
|
||||||
|
|
||||||
!---------------------------------------------------------------------------------------------------
|
|
||||||
!> @brief Calculates IP volume.
|
|
||||||
!> @details The IP volume is calculated differently depending on the cell type.
|
|
||||||
!> 2D cells assume an element depth of one in order to calculate the volume.
|
|
||||||
!> For the hexahedral cell we subdivide the cell into subvolumes of pyramidal
|
|
||||||
!> shape with a cell face as basis and the central ip at the tip. This subvolume is
|
|
||||||
!> calculated as an average of four tetrahedals with three corners on the cell face
|
|
||||||
!> and one corner at the central ip.
|
|
||||||
!---------------------------------------------------------------------------------------------------
|
|
||||||
function IPvolume()
|
|
||||||
|
|
||||||
real(pReal), dimension(theMesh%elem%nIPs,theMesh%nElems) :: IPvolume
|
|
||||||
|
|
||||||
integer :: e,i,c,m,f,n
|
|
||||||
real(pReal), dimension(size(theMesh%elem%cellFace,1),size(theMesh%elem%cellFace,2)) :: subvolume
|
|
||||||
|
|
||||||
c = theMesh%elem%cellType
|
|
||||||
m = size(theMesh%elem%cellFace,1)
|
|
||||||
|
|
||||||
do e = 1,theMesh%nElems
|
|
||||||
select case (c)
|
|
||||||
|
|
||||||
case (1) ! 2D 3node
|
|
||||||
forall (i = 1:theMesh%elem%nIPs) & ! loop over ips=cells in this element
|
|
||||||
IPvolume(i,e) = math_areaTriangle(theMesh%node_0(1:3,mesh_cell2(1,i,e)), &
|
|
||||||
theMesh%node_0(1:3,mesh_cell2(2,i,e)), &
|
|
||||||
theMesh%node_0(1:3,mesh_cell2(3,i,e)))
|
|
||||||
|
|
||||||
case (2) ! 2D 4node
|
|
||||||
forall (i = 1:theMesh%elem%nIPs) & ! loop over ips=cells in this element
|
|
||||||
IPvolume(i,e) = math_areaTriangle(theMesh%node_0(1:3,mesh_cell2(1,i,e)), & ! here we assume a planar shape, so division in two triangles suffices
|
|
||||||
theMesh%node_0(1:3,mesh_cell2(2,i,e)), &
|
|
||||||
theMesh%node_0(1:3,mesh_cell2(3,i,e))) &
|
|
||||||
+ math_areaTriangle(theMesh%node_0(1:3,mesh_cell2(3,i,e)), &
|
|
||||||
theMesh%node_0(1:3,mesh_cell2(4,i,e)), &
|
|
||||||
theMesh%node_0(1:3,mesh_cell2(1,i,e)))
|
|
||||||
|
|
||||||
case (3) ! 3D 4node
|
|
||||||
forall (i = 1:theMesh%elem%nIPs) & ! loop over ips=cells in this element
|
|
||||||
IPvolume(i,e) = math_volTetrahedron(theMesh%node_0(1:3,mesh_cell2(1,i,e)), &
|
|
||||||
theMesh%node_0(1:3,mesh_cell2(2,i,e)), &
|
|
||||||
theMesh%node_0(1:3,mesh_cell2(3,i,e)), &
|
|
||||||
theMesh%node_0(1:3,mesh_cell2(4,i,e)))
|
|
||||||
|
|
||||||
case (4) ! 3D 8node
|
|
||||||
do i = 1,theMesh%elem%nIPs ! loop over ips=cells in this element
|
|
||||||
subvolume = 0.0_pReal
|
|
||||||
forall(f = 1:FE_NipNeighbors(c), n = 1:m) &
|
|
||||||
subvolume(n,f) = math_volTetrahedron(&
|
|
||||||
mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface( n ,f),i,e)), &
|
|
||||||
mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(1+mod(n ,m),f),i,e)), &
|
|
||||||
mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(1+mod(n+1,m),f),i,e)), &
|
|
||||||
mesh_ipCoordinates(1:3,i,e))
|
|
||||||
IPvolume(i,e) = 0.5_pReal * sum(subvolume) ! each subvolume is based on four tetrahedrons, altough the face consists of only two triangles -> averaging factor two
|
|
||||||
enddo
|
|
||||||
|
|
||||||
end select
|
|
||||||
enddo
|
|
||||||
|
|
||||||
end function IPvolume
|
|
||||||
|
|
||||||
|
|
||||||
!---------------------------------------------------------------------------------------------------
|
!---------------------------------------------------------------------------------------------------
|
||||||
!> @brief cell neighborhood
|
!> @brief cell neighborhood
|
||||||
!---------------------------------------------------------------------------------------------------
|
!---------------------------------------------------------------------------------------------------
|
||||||
|
|
@ -1185,74 +1113,6 @@ subroutine mesh_build_ipCoordinates
|
||||||
end subroutine mesh_build_ipCoordinates
|
end subroutine mesh_build_ipCoordinates
|
||||||
|
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
|
||||||
!> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal'
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
|
||||||
subroutine mesh_build_ipAreas(ipArea,ipAreaNormal)
|
|
||||||
|
|
||||||
integer :: e,c,i,f,n,m
|
|
||||||
real(pReal), dimension(theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), intent(out) :: ipArea
|
|
||||||
real(pReal), dimension(3,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), intent(out) :: ipAreaNormal
|
|
||||||
real(pReal), dimension (3,size(theMesh%elem%cellFace,2)) :: nodePos, normals
|
|
||||||
real(pReal), dimension(3) :: normal
|
|
||||||
|
|
||||||
c = theMesh%elem%cellType
|
|
||||||
|
|
||||||
|
|
||||||
do e = 1,theMesh%nElems ! loop over cpElems
|
|
||||||
select case (c)
|
|
||||||
|
|
||||||
case (1,2) ! 2D 3 or 4 node
|
|
||||||
do i = 1,theMesh%elem%nIPs
|
|
||||||
do f = 1,FE_NipNeighbors(c) ! loop over cell faces
|
|
||||||
forall(n = 1: size(theMesh%elem%cellface,1)) &
|
|
||||||
nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(n,f),i,e))
|
|
||||||
normal(1) = nodePos(2,2) - nodePos(2,1) ! x_normal = y_connectingVector
|
|
||||||
normal(2) = -(nodePos(1,2) - nodePos(1,1)) ! y_normal = -x_connectingVector
|
|
||||||
normal(3) = 0.0_pReal
|
|
||||||
ipArea(f,i,e) = norm2(normal)
|
|
||||||
ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
case (3) ! 3D 4node
|
|
||||||
do i = 1,theMesh%elem%nIPs
|
|
||||||
do f = 1,FE_NipNeighbors(c) ! loop over cell faces
|
|
||||||
forall(n = 1: size(theMesh%elem%cellface,1)) &
|
|
||||||
nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(n,f),i,e))
|
|
||||||
normal = math_cross(nodePos(1:3,2) - nodePos(1:3,1), &
|
|
||||||
nodePos(1:3,3) - nodePos(1:3,1))
|
|
||||||
ipArea(f,i,e) = norm2(normal)
|
|
||||||
ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
case (4) ! 3D 8node
|
|
||||||
! for this cell type we get the normal of the quadrilateral face as an average of
|
|
||||||
! four normals of triangular subfaces; since the face consists only of two triangles,
|
|
||||||
! the sum has to be divided by two; this whole prcedure tries to compensate for
|
|
||||||
! probable non-planar cell surfaces
|
|
||||||
m = size(theMesh%elem%cellFace,1)
|
|
||||||
do i = 1,theMesh%elem%nIPs
|
|
||||||
do f = 1,FE_NipNeighbors(c) ! loop over cell faces
|
|
||||||
forall(n = 1: size(theMesh%elem%cellface,1)) &
|
|
||||||
nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(n,f),i,e))
|
|
||||||
forall(n = 1: size(theMesh%elem%cellface,1)) &
|
|
||||||
normals(1:3,n) = 0.5_pReal &
|
|
||||||
* math_cross(nodePos(1:3,1+mod(n ,m)) - nodePos(1:3,n), &
|
|
||||||
nodePos(1:3,1+mod(n+1,m)) - nodePos(1:3,n))
|
|
||||||
normal = 0.5_pReal * sum(normals,2)
|
|
||||||
ipArea(f,i,e) = norm2(normal)
|
|
||||||
ipAreaNormal(1:3,f,i,e) = normal / norm2(normal)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
end select
|
|
||||||
enddo
|
|
||||||
|
|
||||||
end subroutine mesh_build_ipAreas
|
|
||||||
|
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
!> @brief Gives the FE to CP ID mapping by binary search through lookup array
|
!> @brief Gives the FE to CP ID mapping by binary search through lookup array
|
||||||
!! valid questions (what) are 'elem', 'node'
|
!! valid questions (what) are 'elem', 'node'
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue