bugfix: forgot to read file
first draft of nonlocal functions for area normal, area, and volume
This commit is contained in:
Martin Diehl
parent
4e213514bd
commit
dbe15f88f2
|
|
@ -176,7 +176,6 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni
|
|||
microstructureAt, &
|
||||
homogenizationAt
|
||||
|
||||
|
||||
integer :: &
|
||||
fileFormatVersion, &
|
||||
hypoelasticTableStyle, &
|
||||
|
|
@ -194,6 +193,7 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni
|
|||
integer, dimension(:,:), allocatable :: &
|
||||
mapElemSet !< list of elements in elementSet
|
||||
|
||||
inputFile = IO_read_ASCII(trim(modelName)//trim(InputFileExtension))
|
||||
call inputRead_fileFormat(fileFormatVersion, &
|
||||
inputFile)
|
||||
call inputRead_tableStyles(initialcondTableStyle,hypoelasticTableStyle, &
|
||||
|
|
@ -209,10 +209,10 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni
|
|||
call inputRead_mapElemSets(nameElemSet,mapElemSet,&
|
||||
FILEUNIT)
|
||||
|
||||
allocate (mesh_mapFEtoCPelem(2,nElems), source = 0)
|
||||
allocate (mesh_mapFEtoCPelem(2,nElems), source=0)
|
||||
call inputRead_mapElems(hypoelasticTableStyle,nameElemSet,mapElemSet,fileFormatVersion,matNumber,FILEUNIT)
|
||||
|
||||
allocate (mesh_mapFEtoCPnode(2,Nnodes), source = 0)
|
||||
allocate (mesh_mapFEtoCPnode(2,Nnodes), source=0)
|
||||
call inputRead_mapNodes(inputFile)
|
||||
|
||||
call inputRead_elemType(elem, &
|
||||
|
|
@ -423,22 +423,30 @@ subroutine inputRead_mapElems(tableStyle,nameElemSet,mapElemSet,fileFormatVersio
|
|||
character(len=300) :: line, &
|
||||
tmp
|
||||
|
||||
integer, dimension (1+size(mesh_mapFEtoCPelem,2)) :: contInts
|
||||
integer, dimension(:), allocatable :: contInts
|
||||
integer :: i,cpElem
|
||||
|
||||
allocate(contInts(size(mesh_mapFEtoCPelem,2)+1))
|
||||
write(6,*) 'hallo';flush(6)
|
||||
write(6,*) nameElemSet;flush(6)
|
||||
write(6,*) mapElemSet;flush(6)
|
||||
write(6,*) 'map', size(mesh_mapFEtoCPelem,2);flush(6)
|
||||
cpElem = 0
|
||||
contInts = 0
|
||||
rewind(fileUnit)
|
||||
do
|
||||
read (fileUnit,'(A300)',END=620) line
|
||||
write(6,*) trim(line);flush(6)
|
||||
chunkPos = IO_stringPos(line)
|
||||
if (fileFormatVersion < 13) then ! Marc 2016 or earlier
|
||||
if( IO_lc(IO_stringValue(line,chunkPos,1)) == 'hypoelastic' ) then
|
||||
do i=1,3+TableStyle ! skip three (or four if new table style!) lines
|
||||
read (fileUnit,'(A300)') line
|
||||
write(6,*) i;flush(6)
|
||||
enddo
|
||||
contInts = IO_continuousIntValues(fileUnit,size(mesh_mapFEtoCPelem,2),nameElemSet,&
|
||||
mapElemSet,size(nameElemSet))
|
||||
|
||||
exit
|
||||
endif
|
||||
else ! Marc2017 and later
|
||||
|
|
@ -739,7 +747,6 @@ subroutine inputRead_microstructureAndHomogenization(microstructureAt,homogeniza
|
|||
630 end subroutine inputRead_microstructureAndHomogenization
|
||||
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Calculates cell node coordinates from element node coordinates
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -926,6 +933,7 @@ subroutine buildCellNodes(node_cell, &
|
|||
|
||||
end subroutine buildCellNodes
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Calculates IP coordinates as center of cell
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -950,6 +958,135 @@ subroutine buildIPcoordinates(IPcoordinates, &
|
|||
end subroutine buildIPcoordinates
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @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.
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
function IPvolume(elem,node,connectivity)
|
||||
|
||||
type(tElement) :: elem
|
||||
real(pReal), dimension(:,:), intent(in) :: node
|
||||
integer, dimension(:,:,:), intent(in) :: connectivity
|
||||
real(pReal), dimension(elem%nIPs,size(connectivity,3)) :: IPvolume
|
||||
|
||||
integer :: e,i
|
||||
|
||||
do e = 1,size(connectivity,3)
|
||||
do i = 1,elem%nIPs
|
||||
|
||||
select case (elem%cellType)
|
||||
case (1) ! 2D 3node
|
||||
IPvolume(i,e) = math_areaTriangle(node(1:3,connectivity(1,i,e)), &
|
||||
node(1:3,connectivity(2,i,e)), &
|
||||
node(1:3,connectivity(3,i,e)))
|
||||
|
||||
case (2) ! 2D 4node
|
||||
IPvolume(i,e) = math_areaTriangle(node(1:3,connectivity(1,i,e)), & ! assume planar shape, division in two triangles suffices
|
||||
node(1:3,connectivity(2,i,e)), &
|
||||
node(1:3,connectivity(3,i,e))) &
|
||||
+ math_areaTriangle(node(1:3,connectivity(3,i,e)), &
|
||||
node(1:3,connectivity(4,i,e)), &
|
||||
node(1:3,connectivity(1,i,e)))
|
||||
case (3) ! 3D 4node
|
||||
IPvolume(i,e) = math_volTetrahedron(node(1:3,connectivity(1,i,e)), &
|
||||
node(1:3,connectivity(2,i,e)), &
|
||||
node(1:3,connectivity(3,i,e)), &
|
||||
node(1:3,connectivity(4,i,e)))
|
||||
case (4) ! 3D 8node
|
||||
IPvolume(i,e) = math_volTetrahedron(node(1:3,connectivity(1,i,e)), &
|
||||
node(1:3,connectivity(2,i,e)), &
|
||||
node(1:3,connectivity(3,i,e)), &
|
||||
node(1:3,connectivity(6,i,e))) &
|
||||
+ math_volTetrahedron(node(1:3,connectivity(3,i,e)), &
|
||||
node(1:3,connectivity(4,i,e)), &
|
||||
node(1:3,connectivity(1,i,e)), &
|
||||
node(1:3,connectivity(8,i,e))) &
|
||||
+ math_volTetrahedron(node(1:3,connectivity(6,i,e)), &
|
||||
node(1:3,connectivity(7,i,e)), &
|
||||
node(1:3,connectivity(8,i,e)), &
|
||||
node(1:3,connectivity(3,i,e))) &
|
||||
+ math_volTetrahedron(node(1:3,connectivity(8,i,e)), &
|
||||
node(1:3,connectivity(5,i,e)), &
|
||||
node(1:3,connectivity(6,i,e)), &
|
||||
node(1:3,connectivity(1,i,e))) ! first triangulation
|
||||
IPvolume(i,e) = IPvolume(i,e) &
|
||||
+ math_volTetrahedron(node(1:3,connectivity(2,i,e)), &
|
||||
node(1:3,connectivity(3,i,e)), &
|
||||
node(1:3,connectivity(4,i,e)), &
|
||||
node(1:3,connectivity(7,i,e))) &
|
||||
+ math_volTetrahedron(node(1:3,connectivity(4,i,e)), &
|
||||
node(1:3,connectivity(1,i,e)), &
|
||||
node(1:3,connectivity(2,i,e)), &
|
||||
node(1:3,connectivity(5,i,e))) &
|
||||
+ math_volTetrahedron(node(1:3,connectivity(7,i,e)), &
|
||||
node(1:3,connectivity(8,i,e)), &
|
||||
node(1:3,connectivity(5,i,e)), &
|
||||
node(1:3,connectivity(4,i,e))) &
|
||||
+ math_volTetrahedron(node(1:3,connectivity(5,i,e)), &
|
||||
node(1:3,connectivity(6,i,e)), &
|
||||
node(1:3,connectivity(7,i,e)), &
|
||||
node(1:3,connectivity(2,i,e))) ! second triangulation
|
||||
IPvolume(i,e) = IPvolume(i,e) * 0.5_pReal
|
||||
end select
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end function IPvolume
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal'
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine mesh_build_ipAreas(ipArea,ipAreaNormal,elem,nElem,connectivity,node)
|
||||
|
||||
type(tElement) :: elem
|
||||
integer :: nElem
|
||||
integer, dimension(:,:,:), intent(in) :: connectivity
|
||||
real(pReal), dimension(:,:), intent(in) :: node
|
||||
real(pReal), dimension(elem%nIPneighbors,elem%nIPs,nElem), intent(out) :: ipArea
|
||||
real(pReal), dimension(3,elem%nIPneighbors,elem%nIPs,nElem), intent(out) :: ipAreaNormal
|
||||
|
||||
real(pReal), dimension (3,size(elem%cellFace,2)) :: nodePos, normals
|
||||
real(pReal), dimension(3) :: normal
|
||||
integer :: e,i,f,n,m
|
||||
|
||||
m = size(elem%cellFace,1)
|
||||
|
||||
do e = 1,nElem
|
||||
do i = 1,elem%nIPs
|
||||
do f = 1,size(elem%cellFace,2) !????
|
||||
|
||||
forall(n = 1: size(elem%cellface,1)) &
|
||||
nodePos(1:3,n) = node(1:3,connectivity(elem%cellface(n,f),i,e))
|
||||
|
||||
select case (elem%cellType)
|
||||
case (1,2) ! 2D 3 or 4 node
|
||||
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
|
||||
case (3) ! 3D 4node
|
||||
normal = math_cross(nodePos(1:3,2) - nodePos(1:3,1), &
|
||||
nodePos(1:3,3) - nodePos(1:3,1))
|
||||
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
|
||||
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)
|
||||
end select
|
||||
|
||||
ipArea(f,i,e) = norm2(normal)
|
||||
ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end subroutine mesh_build_ipAreas
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Gives the FE to CP ID mapping by binary search through lookup array
|
||||
!! valid questions (what) are 'elem', 'node'
|
||||
|
|
|
|||
Loading…
Reference in New Issue