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

major debugging and streamlining of the calculation of ipVolumes and 

ipAreas
This commit is contained in:
Philip Eisenlohr 2009-01-19 18:42:31 +00:00
parent e3e42fdf44
commit 82e65648a2
1 changed files with 153 additions and 53 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

206
trunk/mesh.f90
View File

@ -51,6 +51,7 @@
integer(pInt), dimension(:,:), allocatable, target :: mesh_mapFEtoCPelem,mesh_mapFEtoCPnode integer(pInt), dimension(:,:), allocatable, target :: mesh_mapFEtoCPelem,mesh_mapFEtoCPnode
integer(pInt), dimension(:,:), allocatable :: mesh_element, mesh_sharedElem integer(pInt), dimension(:,:), allocatable :: mesh_element, mesh_sharedElem
integer(pInt), dimension(:,:,:,:), allocatable :: mesh_ipNeighborhood integer(pInt), dimension(:,:,:,:), allocatable :: mesh_ipNeighborhood
real(pReal), dimension(:,:,:), allocatable :: mesh_subNodeCoord ! coordinates of subnodes per element real(pReal), dimension(:,:,:), allocatable :: mesh_subNodeCoord ! coordinates of subnodes per element
real(pReal), dimension(:,:), allocatable :: mesh_ipVolume ! volume associated with IP real(pReal), dimension(:,:), allocatable :: mesh_ipVolume ! volume associated with IP
real(pReal), dimension(:,:,:), allocatable :: mesh_ipArea ! area of interface to neighboring IP real(pReal), dimension(:,:,:), allocatable :: mesh_ipArea ! area of interface to neighboring IP
@ -91,11 +92,11 @@
/) /)
integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_NsubNodes = & integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_NsubNodes = &
(/19, & ! element 7 (/19, & ! element 7
0, & ! element 134 0, & ! element 134
0, & ! element 11 0, & ! element 11
0, & ! element 27 0, & ! element 27
0, & ! element 157 0, & ! element 157
0 & ! element 136 0 & ! element 136
/) /)
integer(pInt), dimension(FE_maxNipNeighbors,FE_Nelemtypes), parameter :: FE_NfaceNodes = & integer(pInt), dimension(FE_maxNipNeighbors,FE_Nelemtypes), parameter :: FE_NfaceNodes = &
reshape((/& reshape((/&
@ -225,7 +226,7 @@
1, 2, 0, 0, 0, 0, 0, 0, & ! element 7 1, 2, 0, 0, 0, 0, 0, 0, & ! element 7
2, 3, 0, 0, 0, 0, 0, 0, & 2, 3, 0, 0, 0, 0, 0, 0, &
3, 4, 0, 0, 0, 0, 0, 0, & 3, 4, 0, 0, 0, 0, 0, 0, &
4, 5, 0, 0, 0, 0, 0, 0, & 4, 1, 0, 0, 0, 0, 0, 0, &
1, 5, 0, 0, 0, 0, 0, 0, & 1, 5, 0, 0, 0, 0, 0, 0, &
2, 6, 0, 0, 0, 0, 0, 0, & 2, 6, 0, 0, 0, 0, 0, 0, &
3, 7, 0, 0, 0, 0, 0, 0, & 3, 7, 0, 0, 0, 0, 0, 0, &
@ -343,7 +344,7 @@
1,13,25,12, & 1,13,25,12, &
12,25,27,21, & 12,25,27,21, &
1, 9,22,13, & 1, 9,22,13, &
13,22,27,26, & 13,22,27,25, &
1,12,21, 9, & 1,12,21, 9, &
2,10,23,14, & ! 2,10,23,14, & !
9,22,27,21, & 9,22,27,21, &
@ -370,7 +371,7 @@
5,17,26,20, & 5,17,26,20, &
13,25,27,22, & 13,25,27,22, &
6,14,23,18, & ! 6,14,23,18, & !
17,26,27,18, & 17,26,27,22, &
18,23,27,26, & 18,23,27,26, &
6,17,22,14, & 6,17,22,14, &
6,18,26,17, & 6,18,26,17, &
@ -393,10 +394,10 @@
0, 0, 0, 0, & 0, 0, 0, 0, &
0, 0, 0, 0, & 0, 0, 0, 0, &
0, 0, 0, 0, & 0, 0, 0, 0, &
1, 3, 2, 0, & ! element 134 1, 1, 3, 2, & ! element 134
1, 2, 4, 0, & 1, 1, 2, 4, &
2, 3, 4, 0, & 2, 2, 3, 4, &
1, 4, 3, 0, & 1, 1, 4, 3, &
0, 0, 0, 0, & 0, 0, 0, 0, &
0, 0, 0, 0, & 0, 0, 0, 0, &
0, 0, 0, 0, & ! 0, 0, 0, 0, & !
@ -672,16 +673,20 @@
! function mesh_build_ipNeighorhood() ! function mesh_build_ipNeighorhood()
! --------------------------- ! ---------------------------
!*********************************************************** !***********************************************************
! initialization ! initialization
!*********************************************************** !***********************************************************
SUBROUTINE mesh_init () SUBROUTINE mesh_init ()
use prec, only: pInt use prec, only: pInt
use IO, only: IO_error,IO_open_InputFile use IO, only: IO_error,IO_open_InputFile
use FEsolving, only: FE_get_solverSymmetry use FEsolving, only: FE_get_solverSymmetry
implicit none implicit none
integer(pInt), parameter :: fileUnit = 222 integer(pInt), parameter :: fileUnit = 222
mesh_Nelems = 0_pInt mesh_Nelems = 0_pInt
mesh_NcpElems = 0_pInt mesh_NcpElems = 0_pInt
mesh_Nnodes = 0_pInt mesh_Nnodes = 0_pInt
@ -693,8 +698,11 @@
mesh_NelemSets = 0_pInt mesh_NelemSets = 0_pInt
mesh_maxNelemInSet = 0_pInt mesh_maxNelemInSet = 0_pInt
! call to various subroutines to parse the stuff from the input file... ! call to various subroutines to parse the stuff from the input file...
if (IO_open_inputFile(fileUnit)) then if (IO_open_inputFile(fileUnit)) then
call FE_get_solverSymmetry(fileUnit) call FE_get_solverSymmetry(fileUnit)
call mesh_get_meshDimensions(fileUnit) call mesh_get_meshDimensions(fileUnit)
call mesh_build_nodeMapping(fileUnit) call mesh_build_nodeMapping(fileUnit)
@ -717,10 +725,12 @@
END SUBROUTINE END SUBROUTINE
!*********************************************************** !***********************************************************
! mapping of FE element types to internal representation ! mapping of FE element types to internal representation
!*********************************************************** !***********************************************************
FUNCTION FE_mapElemtype(what) FUNCTION FE_mapElemtype(what)
implicit none implicit none
character(len=*), intent(in) :: what character(len=*), intent(in) :: what
@ -742,8 +752,11 @@
case default case default
FE_mapElemtype = 0 ! unknown element --> should raise an error upstream..! FE_mapElemtype = 0 ! unknown element --> should raise an error upstream..!
end select end select
END FUNCTION END FUNCTION
!*********************************************************** !***********************************************************
! FE to CP id mapping by binary search thru lookup array ! FE to CP id mapping by binary search thru lookup array
! !
@ -798,12 +811,15 @@
END FUNCTION END FUNCTION
!*********************************************************** !***********************************************************
! find face-matching element of same type ! find face-matching element of same type
!!*********************************************************** !!***********************************************************
FUNCTION mesh_faceMatch(face,elem) FUNCTION mesh_faceMatch(face,elem)
use prec, only: pInt use prec, only: pInt
implicit none implicit none
integer(pInt) face,elem integer(pInt) face,elem
integer(pInt) mesh_faceMatch integer(pInt) mesh_faceMatch
integer(pInt), dimension(FE_NfaceNodes(face,mesh_element(2,elem))) :: nodeMap integer(pInt), dimension(FE_NfaceNodes(face,mesh_element(2,elem))) :: nodeMap
@ -812,6 +828,7 @@
minN = mesh_maxNsharedElems+1 ! init to worst case minN = mesh_maxNsharedElems+1 ! init to worst case
mesh_faceMatch = 0_pInt ! intialize to "no match found" mesh_faceMatch = 0_pInt ! intialize to "no match found"
t = mesh_element(2,elem) ! figure elemType t = mesh_element(2,elem) ! figure elemType
do faceNode=1,FE_NfaceNodes(face,t) ! loop over nodes on face do faceNode=1,FE_NfaceNodes(face,t) ! loop over nodes on face
nodeMap(faceNode) = mesh_FEasCP('node',mesh_element(4+FE_nodeOnFace(faceNode,face,t),elem)) ! CP id of face node nodeMap(faceNode) = mesh_FEasCP('node',mesh_element(4+FE_nodeOnFace(faceNode,face,t),elem)) ! CP id of face node
NsharedElems = mesh_sharedElem(1,nodeMap(faceNode)) ! figure # shared elements for this node NsharedElems = mesh_sharedElem(1,nodeMap(faceNode)) ! figure # shared elements for this node
@ -841,6 +858,7 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
END FUNCTION END FUNCTION
!******************************************************************** !********************************************************************
! get count of elements, nodes, and cp elements in mesh ! get count of elements, nodes, and cp elements in mesh
! for subsequent array allocations ! for subsequent array allocations
@ -849,16 +867,21 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
! _Nelems, _Nnodes, _NcpElems ! _Nelems, _Nnodes, _NcpElems
!******************************************************************** !********************************************************************
SUBROUTINE mesh_get_meshDimensions (unit) SUBROUTINE mesh_get_meshDimensions (unit)
use prec, only: pInt use prec, only: pInt
use IO use IO
implicit none implicit none
integer(pInt) unit,i,pos(41) integer(pInt) unit,i,pos(41)
character*300 line character*300 line
610 FORMAT(A300) 610 FORMAT(A300)
rewind(unit) rewind(unit)
do do
read (unit,610,END=620) line read (unit,610,END=620) line
pos = IO_stringPos(line,20) pos = IO_stringPos(line,20)
select case ( IO_lc(IO_StringValue(line,pos,1))) select case ( IO_lc(IO_StringValue(line,pos,1)))
case('table') case('table')
if (pos(1) == 6) then if (pos(1) == 6) then
@ -880,10 +903,13 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
end do end do
mesh_NcpElems = mesh_NcpElems + IO_countContinousIntValues(unit) mesh_NcpElems = mesh_NcpElems + IO_countContinousIntValues(unit)
end select end select
end do end do
620 return 620 return
END SUBROUTINE END SUBROUTINE
!!******************************************************************** !!********************************************************************
! get maximum count of nodes, IPs, IP neighbors, and shared elements ! get maximum count of nodes, IPs, IP neighbors, and shared elements
@ -924,7 +950,9 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(t)) mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(t))
mesh_maxNips = max(mesh_maxNips,FE_Nips(t)) mesh_maxNips = max(mesh_maxNips,FE_Nips(t))
mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(t)) mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(t))
mesh_maxNsubNodes = max(mesh_maxNsubNodes,FE_NsubNodes(t)) mesh_maxNsubNodes = max(mesh_maxNsubNodes,FE_NsubNodes(t))
node_seen = 0_pInt node_seen = 0_pInt
do j=1,FE_Nnodes(t) do j=1,FE_Nnodes(t)
n = mesh_FEasCP('node',IO_IntValue (line,pos,j+2)) n = mesh_FEasCP('node',IO_IntValue (line,pos,j+2))
@ -950,16 +978,22 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
! allocate globals: mesh_nameElemSet, mesh_mapElemSet ! allocate globals: mesh_nameElemSet, mesh_mapElemSet
!******************************************************************** !********************************************************************
SUBROUTINE mesh_build_elemSetMapping (unit) SUBROUTINE mesh_build_elemSetMapping (unit)
use prec, only: pInt use prec, only: pInt
use IO use IO
implicit none implicit none
integer unit, elem_set integer unit, elem_set
character*300 line character*300 line
integer(pInt), dimension (9) :: pos ! count plus 4 entities on a line integer(pInt), dimension (9) :: pos ! count plus 4 entities on a line
610 FORMAT(A300) 610 FORMAT(A300)
allocate (mesh_nameElemSet(mesh_NelemSets)) allocate (mesh_nameElemSet(mesh_NelemSets))
allocate (mesh_mapElemSet(1+mesh_maxNelemInSet,mesh_NelemSets)) ; mesh_mapElemSet = 0_pInt allocate (mesh_mapElemSet(1+mesh_maxNelemInSet,mesh_NelemSets)) ; mesh_mapElemSet = 0_pInt
elem_set = 0_pInt elem_set = 0_pInt
rewind(unit) rewind(unit)
do do
read (unit,610,END=620) line read (unit,610,END=620) line
@ -971,8 +1005,10 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
mesh_mapElemSet(:,elem_set) = IO_continousIntValues(unit,mesh_maxNelemInSet,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) mesh_mapElemSet(:,elem_set) = IO_continousIntValues(unit,mesh_maxNelemInSet,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets)
end if end if
end do end do
620 return 620 return
END SUBROUTINE END SUBROUTINE
!******************************************************************** !********************************************************************
! Build node mapping from FEM to CP ! Build node mapping from FEM to CP
@ -981,17 +1017,22 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
! _mapFEtoCPnode ! _mapFEtoCPnode
!******************************************************************** !********************************************************************
SUBROUTINE mesh_build_nodeMapping (unit) SUBROUTINE mesh_build_nodeMapping (unit)
use prec, only: pInt use prec, only: pInt
use math, only: qsort use math, only: qsort
use IO use IO
implicit none implicit none
integer(pInt), dimension (mesh_Nnodes) :: node_count integer(pInt), dimension (mesh_Nnodes) :: node_count
integer(pInt) unit,i integer(pInt) unit,i
integer(pInt), dimension (133) :: pos integer(pInt), dimension (133) :: pos
character*300 line character*300 line
610 FORMAT(A300) 610 FORMAT(A300)
allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt
node_count(:) = 0_pInt node_count(:) = 0_pInt
rewind(unit) rewind(unit)
do do
read (unit,610,END=620) line read (unit,610,END=620) line
@ -1006,10 +1047,13 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
exit exit
end if end if
end do end do
620 call qsort(mesh_mapFEtoCPnode,1,size(mesh_mapFEtoCPnode,2)) 620 call qsort(mesh_mapFEtoCPnode,1,size(mesh_mapFEtoCPnode,2))
return return
END SUBROUTINE END SUBROUTINE
!******************************************************************** !********************************************************************
! Build element mapping from FEM to CP ! Build element mapping from FEM to CP
! !
@ -1017,18 +1061,24 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
! _mapFEtoCPelem ! _mapFEtoCPelem
!******************************************************************** !********************************************************************
SUBROUTINE mesh_build_elemMapping (unit) SUBROUTINE mesh_build_elemMapping (unit)
use prec, only: pInt use prec, only: pInt
use math, only: qsort use math, only: qsort
use IO use IO
implicit none implicit none
integer unit, i,CP_elem integer unit, i,CP_elem
character*300 line character*300 line
integer(pInt), dimension (3) :: pos integer(pInt), dimension (3) :: pos
integer(pInt), dimension (1+mesh_NcpElems) :: contInts integer(pInt), dimension (1+mesh_NcpElems) :: contInts
610 FORMAT(A300) 610 FORMAT(A300)
allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt
CP_elem = 0_pInt CP_elem = 0_pInt
rewind(unit) rewind(unit)
do do
read (unit,610,END=620) line read (unit,610,END=620) line
@ -1045,10 +1095,13 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
enddo enddo
end if end if
end do end do
620 call qsort(mesh_mapFEtoCPelem,1,size(mesh_mapFEtoCPelem,2)) ! should be mesh_NcpElems 620 call qsort(mesh_mapFEtoCPelem,1,size(mesh_mapFEtoCPelem,2)) ! should be mesh_NcpElems
return return
END SUBROUTINE END SUBROUTINE
!******************************************************************** !********************************************************************
! store x,y,z coordinates of all nodes in mesh ! store x,y,z coordinates of all nodes in mesh
! !
@ -1056,16 +1109,21 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
! _node ! _node
!******************************************************************** !********************************************************************
SUBROUTINE mesh_build_nodes (unit) SUBROUTINE mesh_build_nodes (unit)
use prec, only: pInt use prec, only: pInt
use IO use IO
implicit none implicit none
integer unit,i,j,m integer unit,i,j,m
integer(pInt), dimension(3) :: pos integer(pInt), dimension(3) :: pos
integer(pInt), dimension(5), parameter :: node_ends = (/0,10,30,50,70/) integer(pInt), dimension(5), parameter :: node_ends = (/0,10,30,50,70/)
character*300 line character*300 line
allocate ( mesh_node (3,mesh_Nnodes) ) allocate ( mesh_node (3,mesh_Nnodes) )
mesh_node(:,:) = 0_pInt mesh_node(:,:) = 0_pInt
610 FORMAT(A300) 610 FORMAT(A300)
rewind(unit) rewind(unit)
do do
read (unit,610,END=620) line read (unit,610,END=620) line
@ -1082,9 +1140,12 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
exit exit
end if end if
end do end do
620 return 620 return
END SUBROUTINE END SUBROUTINE
!******************************************************************** !********************************************************************
! store FEid, type, mat, tex, and node list per element ! store FEid, type, mat, tex, and node list per element
! !
@ -1092,16 +1153,21 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
! _element ! _element
!******************************************************************** !********************************************************************
SUBROUTINE mesh_build_elements (unit) SUBROUTINE mesh_build_elements (unit)
use prec, only: pInt use prec, only: pInt
use IO use IO
implicit none implicit none
integer unit,i,j,sv,val,CP_elem integer unit,i,j,sv,val,CP_elem
integer(pInt), dimension(133) :: pos integer(pInt), dimension(133) :: pos
integer(pInt), dimension(1+mesh_NcpElems) :: contInts integer(pInt), dimension(1+mesh_NcpElems) :: contInts
character*300 line character*300 line
allocate (mesh_element (4+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt allocate (mesh_element (4+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt
610 FORMAT(A300) 610 FORMAT(A300)
rewind(unit) rewind(unit)
do do
read (unit,610,END=620) line read (unit,610,END=620) line
@ -1158,8 +1224,11 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
read (unit,610,END=620) line read (unit,610,END=620) line
endif endif
enddo enddo
620 return 620 return
END SUBROUTINE END SUBROUTINE
!******************************************************************** !********************************************************************
! build list of elements shared by each node in mesh ! build list of elements shared by each node in mesh
@ -1168,17 +1237,22 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
! _sharedElem ! _sharedElem
!******************************************************************** !********************************************************************
SUBROUTINE mesh_build_sharedElems (unit) SUBROUTINE mesh_build_sharedElems (unit)
use prec, only: pInt use prec, only: pInt
use IO use IO
implicit none implicit none
integer(pint) unit,i,j,n,e integer(pint) unit,i,j,n,e
integer(pInt), dimension (133) :: pos integer(pInt), dimension (133) :: pos
integer(pInt), dimension (:), allocatable :: node_seen integer(pInt), dimension (:), allocatable :: node_seen
character*300 line character*300 line
610 FORMAT(A300) 610 FORMAT(A300)
allocate(node_seen(maxval(FE_Nnodes))) allocate(node_seen(maxval(FE_Nnodes)))
allocate ( mesh_sharedElem( 1+mesh_maxNsharedElems,mesh_Nnodes) ) allocate ( mesh_sharedElem( 1+mesh_maxNsharedElems,mesh_Nnodes) )
mesh_sharedElem(:,:) = 0_pInt mesh_sharedElem(:,:) = 0_pInt
rewind(unit) rewind(unit)
do do
read (unit,610,END=620) line read (unit,610,END=620) line
@ -1204,9 +1278,12 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
exit exit
end if end if
end do end do
620 return 620 return
END SUBROUTINE END SUBROUTINE
!*********************************************************** !***********************************************************
! build up of IP neighborhood ! build up of IP neighborhood
! !
@ -1220,6 +1297,7 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
integer(pInt) e,t,i,j,k,n integer(pInt) e,t,i,j,k,n
integer(pInt) neighbor,neighboringElem,neighboringIP,matchingElem,faceNode,linkingNode integer(pInt) neighbor,neighboringElem,neighboringIP,matchingElem,faceNode,linkingNode
allocate(mesh_ipNeighborhood(2,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipNeighborhood = 0_pInt allocate(mesh_ipNeighborhood(2,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipNeighborhood = 0_pInt
do e = 1,mesh_NcpElems ! loop over cpElems do e = 1,mesh_NcpElems ! loop over cpElems
@ -1258,9 +1336,10 @@ matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matc
end do end do
return return
END SUBROUTINE END SUBROUTINE
!*********************************************************** !***********************************************************
! assignment of coordinates for subnodes in each cp element ! assignment of coordinates for subnodes in each cp element
@ -1282,15 +1361,15 @@ matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matc
do n = 1,FE_Nnodes(t) do n = 1,FE_Nnodes(t)
mesh_subNodeCoord(:,n,e) = mesh_node(:,mesh_FEasCP('node',mesh_element(4+n,e))) ! loop over nodes of this element type mesh_subNodeCoord(:,n,e) = mesh_node(:,mesh_FEasCP('node',mesh_element(4+n,e))) ! loop over nodes of this element type
enddo enddo
do n = 1,FE_NsubNodes(t) ! now for the treu subnodes do n = 1,FE_NsubNodes(t) ! now for the true subnodes
do p = 1,mesh_maxNnodes ! loop through parents do p = 1,FE_Nnodes(t) ! loop through parents
if (FE_subNodeParent(p,n,t) > 0) & ! valid parent node if (FE_subNodeParent(p,n,t) > 0) & ! valid parent node
mesh_subNodeCoord(:,n+FE_Nnodes(t),e) = & mesh_subNodeCoord(:,n+FE_Nnodes(t),e) = &
mesh_subNodeCoord(:,n+FE_Nnodes(t),e) + & mesh_subNodeCoord(:,n+FE_Nnodes(t),e) + &
mesh_node(:,mesh_FEasCP('node',mesh_element(4+FE_subNodeParent(p,n,t),e))) ! add up parents mesh_node(:,mesh_FEasCP('node',mesh_element(4+FE_subNodeParent(p,n,t),e))) ! add up parents
enddo enddo
mesh_subNodeCoord(:,n+FE_Nnodes(t),e) = mesh_subNodeCoord(:,n+FE_Nnodes(t),e) / count(FE_subNodeParent(:,n,t) > 0) mesh_subNodeCoord(:,n+FE_Nnodes(t),e) = mesh_subNodeCoord(:,n+FE_Nnodes(t),e) / count(FE_subNodeParent(:,n,t) > 0)
enddo enddo
enddo enddo
return return
@ -1311,9 +1390,11 @@ matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matc
implicit none implicit none
integer(pInt) e,f,t,i,j,k,n integer(pInt) e,f,t,i,j,k,n
integer(pInt), parameter :: Ntriangles = FE_NipFaceNodes-2 ! each interface is made up of this many triangles
integer(pInt), dimension(mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNode ! flagList to find subnodes determining center of grav integer(pInt), dimension(mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNode ! flagList to find subnodes determining center of grav
real(pReal), dimension(3,mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNodePos ! coordinates of subnodes determining center of grav real(pReal), dimension(3,mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNodePos ! coordinates of subnodes determining center of grav
real(pReal), dimension (3,FE_NipFaceNodes) :: nPos ! coordinates of nodes on IP face real(pReal), dimension (3,FE_NipFaceNodes) :: nPos ! coordinates of nodes on IP face
real(pReal), dimension(Ntriangles,FE_NipFaceNodes) :: volume ! volumes of possible tetrahedra
real(pReal), dimension(3) :: centerOfGravity real(pReal), dimension(3) :: centerOfGravity
allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems)) ; mesh_ipVolume = 0.0_pReal allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems)) ; mesh_ipVolume = 0.0_pReal
@ -1323,35 +1404,36 @@ matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matc
do i = 1,FE_Nips(t) ! loop over IPs of elem do i = 1,FE_Nips(t) ! loop over IPs of elem
gravityNode = 0_pInt ! reset flagList gravityNode = 0_pInt ! reset flagList
gravityNodePos = 0.0_pReal ! reset coordinates gravityNodePos = 0.0_pReal ! reset coordinates
do f = 1,FE_NipNeighbors(t) ! loop over neighbors of IP do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP
do n = 1,FE_NipFaceNodes ! loop over nodes on interface do n = 1,FE_NipFaceNodes ! loop over nodes on interface
gravityNode(FE_subNodeOnIPFace(n,f,i,t)) = 1 gravityNode(FE_subNodeOnIPFace(n,f,i,t)) = 1
gravityNodePos(:,FE_subNodeOnIPFace(n,f,i,t)) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) gravityNodePos(:,FE_subNodeOnIPFace(n,f,i,t)) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e)
enddo enddo
enddo enddo
do j = 1,8+FE_maxNsubNodes-1 ! walk through entire flagList except last do j = 1,mesh_maxNnodes+mesh_maxNsubNodes-1 ! walk through entire flagList except last
if (gravityNode(j) > 0_pInt) then ! valid node index if (gravityNode(j) > 0_pInt) then ! valid node index
do k = j,8+FE_maxNsubNodes ! walk through remainder of list do k = j+1,mesh_maxNnodes+mesh_maxNsubNodes ! walk through remainder of list
if (all((mesh_subNodeCoord(:,gravityNode(j),e) - & if (all((gravityNodePos(:,j) - gravityNodePos(:,k)) == 0.0_pReal)) then ! found match
mesh_subNodeCoord(:,gravityNode(k),e)) == 0.0_pReal)) then ! found match gravityNode(j) = 0_pInt ! delete first instance
gravityNode(j) = 0_pInt ! delete first instance
gravityNodePos(:,j) = 0.0_pReal gravityNodePos(:,j) = 0.0_pReal
exit ! continue with next suspect exit ! continue with next suspect
endif endif
enddo enddo
endif endif
enddo enddo
centerOfGravity = sum(gravityNodePos,2)/count(gravityNode > 0) centerOfGravity = sum(gravityNodePos,2)/count(gravityNode > 0)
do f = 1,FE_NipNeighbors(t) ! loop over neighbors of IP and add tetrahedra which connect to CoG do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP and add tetrahedra which connect to CoG
forall (n=1:FE_NipFaceNodes) nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) forall (n = 1:FE_NipFaceNodes) nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e)
mesh_ipVolume(i,e) = mesh_ipVolume(i,e) + math_volTetrahedron(nPos(:,1),nPos(:,2),nPos(:,3),centerOfGravity) forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles) & ! start at each interface node and build valid triangles to cover interface
mesh_ipVolume(i,e) = mesh_ipVolume(i,e) + math_volTetrahedron(nPos(:,1),nPos(:,3),nPos(:,4),centerOfGravity) volume(j,n) = math_volTetrahedron(nPos(:,n), & ! calc volume of respective tetrahedron to CoG
mesh_ipVolume(i,e) = mesh_ipVolume(i,e) + math_volTetrahedron(nPos(:,2),nPos(:,3),nPos(:,4),centerOfGravity) nPos(:,1+mod(n+j-1,FE_NipFaceNodes)), &
mesh_ipVolume(i,e) = mesh_ipVolume(i,e) + math_volTetrahedron(nPos(:,2),nPos(:,4),nPos(:,1),centerOfGravity) nPos(:,1+mod(n+j-0,FE_NipFaceNodes)), &
centerOfGravity)
mesh_ipVolume(i,e) = mesh_ipVolume(i,e) + sum(volume) ! add contribution from this interface
enddo enddo
mesh_ipVolume(i,e) = 0.5_pReal * mesh_ipVolume(i,e) mesh_ipVolume(i,e) = mesh_ipVolume(i,e) / FE_NipFaceNodes ! renormalize with interfaceNodeNum due to loop over them
enddo enddo
enddo enddo
return return
@ -1372,29 +1454,29 @@ matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matc
implicit none implicit none
integer(pInt) e,f,t,i,j,k,n integer(pInt) e,f,t,i,j,k,n
real(pReal), dimension (3,FE_NipFaceNodes) :: nPos ! coordinates of nodes on IP face integer(pInt), parameter :: Ntriangles = FE_NipFaceNodes-2 ! each interface is made up of this many triangles
real(pReal), dimension(3,4) :: normal real(pReal), dimension (3,FE_NipFaceNodes) :: nPos ! coordinates of nodes on IP face
real(pReal), dimension(4) :: area real(pReal), dimension(3,Ntriangles,FE_NipFaceNodes) :: normal
real(pReal), dimension(Ntriangles,FE_NipFaceNodes) :: area
allocate(mesh_ipArea(mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipArea = 0.0_pReal allocate(mesh_ipArea(mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipArea = 0.0_pReal
allocate(mesh_ipAreaNormal(3,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipAreaNormal = 0.0_pReal allocate(mesh_ipAreaNormal(3,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipAreaNormal = 0.0_pReal
do e = 1,mesh_NcpElems ! loop over cpElems do e = 1,mesh_NcpElems ! loop over cpElems
t = mesh_element(2,e) ! get elemType t = mesh_element(2,e) ! get elemType
do i = 1,FE_Nips(t) ! loop over IPs of elem do i = 1,FE_Nips(t) ! loop over IPs of elem
do f = 1,FE_NipNeighbors(t) ! loop over neighbors of IP do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP
forall (n = 1:FE_NipFaceNodes) & ! loop over nodes on interface forall (n = 1:FE_NipFaceNodes) nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e)
nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) ! get shorthand name forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles) ! start at each interface node and build valid triangles to cover interface
normal(:,1) = math_vectorproduct(nPos(:,2)-nPos(:,1),nPos(:,3)-nPos(:,1)) normal(:,j,n) = math_vectorproduct(nPos(:,1+mod(n+j-1,FE_NipFaceNodes)) - nPos(:,n), & ! calc their normal vectors
normal(:,2) = math_vectorproduct(nPos(:,3)-nPos(:,1),nPos(:,4)-nPos(:,1)) nPos(:,1+mod(n+j-0,FE_NipFaceNodes)) - nPos(:,n))
normal(:,3) = math_vectorproduct(nPos(:,3)-nPos(:,2),nPos(:,4)-nPos(:,2)) area(j,n) = dsqrt(sum(normal(:,j,n)*normal(:,j,n))) ! and area
normal(:,4) = math_vectorproduct(nPos(:,4)-nPos(:,2),nPos(:,1)-nPos(:,2)) end forall
forall (n = 1:4) forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles, area(j,n) > 0.0_pReal) &
area(n) = dsqrt(sum(normal(:,n)*normal(:,n))) ! length of normal vectors (i.e. area) normal(:,j,n) = normal(:,j,n) / area(j,n) ! make unit normal
normal(:,n) = normal(:,n) / area(n) ! make unit normal
endforall mesh_ipArea(f,i,e) = sum(area) / (FE_NipFaceNodes*2.0_pReal) ! area of parallelograms instead of triangles
mesh_ipArea(f,i,e) = sum(area) / 2.0_pReal ! counted area twice mesh_ipAreaNormal(:,f,i,e) = sum(sum(normal,3),2) / count(area > 0.0_pReal) ! average of all valid normals
mesh_ipAreaNormal(:,f,i,e) = sum(normal,2) / 4.0_pReal ! average of all four possible normals
enddo enddo
enddo enddo
enddo enddo
@ -1403,20 +1485,22 @@ matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matc
END SUBROUTINE END SUBROUTINE
!*********************************************************** !***********************************************************
! write statistics regarding input file parsing ! write statistics regarding input file parsing
! to the output file ! to the output file
! !
!*********************************************************** !***********************************************************
SUBROUTINE mesh_tell_statistics() SUBROUTINE mesh_tell_statistics()
use prec, only: pInt use prec, only: pInt
use IO, only: IO_error use IO, only: IO_error
implicit none implicit none
integer(pInt), dimension (:,:), allocatable :: mesh_MatTex integer(pInt), dimension (:,:), allocatable :: mesh_MatTex
character(len=64) fmt character(len=64) fmt
integer(pInt) i
integer(pInt) i,e,f
if (mesh_maxValStateVar(1) == 0) call IO_error(110) ! no materials specified if (mesh_maxValStateVar(1) == 0) call IO_error(110) ! no materials specified
if (mesh_maxValStateVar(2) == 0) call IO_error(120) ! no textures specified if (mesh_maxValStateVar(2) == 0) call IO_error(120) ! no textures specified
@ -1429,6 +1513,21 @@ matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matc
enddo enddo
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write (6,*)
write (6,*) "Input Parser: IP NEIGHBORHOOD"
write (6,*)
write (6,"(a13,x,e15.8)") "total volume", sum(mesh_ipVolume)
write (6,*)
write (6,"(a5,x,a5,x,a15,x,a5,x,a15,x,a16)") "elem","IP","volume","face","area","-- normal --"
do e = 1,mesh_NcpElems
do i = 1,FE_Nips(mesh_element(2,e))
write (6,"(i5,x,i5,x,e15.8)") e,i,mesh_IPvolume(i,e)
do f = 1,FE_NipNeighbors(mesh_element(2,e))
! write (6,"(i33,x,e15.8,x,3(f6.3,x))") f,mesh_ipArea(f,i,e),mesh_ipAreaNormal(:,f,i,e)
enddo
enddo
enddo
write (6,*) write (6,*)
write (6,*) "Input Parser: STATISTICS" write (6,*) "Input Parser: STATISTICS"
write (6,*) write (6,*)
@ -1453,6 +1552,7 @@ matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matc
write (6,*) write (6,*)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
return return
END SUBROUTINE END SUBROUTINE