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

reworked build_ipNeighborhood 

added C3D8R reduced integration hexahedral element
This commit is contained in:
Philip Eisenlohr 2010-05-10 14:54:59 +00:00
parent 2405a51042
commit 0416c9a616
1 changed files with 243 additions and 190 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

433
code/mesh.f90
View File

@ -35,6 +35,7 @@
! _NfaceNodes : # nodes per face in a specific type of element ! _NfaceNodes : # nodes per face in a specific type of element
! _nodeOnFace : list of node indices on each face of a specific type of element ! _nodeOnFace : list of node indices on each face of a specific type of element
! _maxNnodesAtIP : max number of (equivalent) nodes attached to an IP
! _nodesAtIP : map IP index to two node indices in a specific type of element ! _nodesAtIP : map IP index to two node indices in a specific type of element
! _ipNeighborhood : 6 or less neighboring IPs as [element_num, IP_index] ! _ipNeighborhood : 6 or less neighboring IPs as [element_num, IP_index]
! _NsubNodes : # subnodes required to fully define all IP volumes ! _NsubNodes : # subnodes required to fully define all IP volumes
@ -60,18 +61,19 @@
real(pReal), dimension(:,:,:,:), allocatable :: mesh_ipAreaNormal ! area normal of interface to neighboring IP real(pReal), dimension(:,:,:,:), allocatable :: mesh_ipAreaNormal ! area normal of interface to neighboring IP
real(pReal), allocatable :: mesh_node (:,:) real(pReal), allocatable :: mesh_node (:,:)
integer(pInt), dimension(:,:,:,:), allocatable :: FE_nodesAtIP integer(pInt), dimension(:,:,:), allocatable :: FE_nodesAtIP
integer(pInt), dimension(:,:,:), allocatable :: FE_ipNeighbor integer(pInt), dimension(:,:,:), allocatable :: FE_ipNeighbor
integer(pInt), dimension(:,:,:), allocatable :: FE_subNodeParent integer(pInt), dimension(:,:,:), allocatable :: FE_subNodeParent
integer(pInt), dimension(:,:,:,:), allocatable :: FE_subNodeOnIPFace integer(pInt), dimension(:,:,:,:), allocatable :: FE_subNodeOnIPFace
integer(pInt) :: hypoelasticTableStyle integer(pInt) :: hypoelasticTableStyle
integer(pInt) :: initialcondTableStyle integer(pInt) :: initialcondTableStyle
integer(pInt), parameter :: FE_Nelemtypes = 7 integer(pInt), parameter :: FE_Nelemtypes = 8
integer(pInt), parameter :: FE_maxNnodes = 8 integer(pInt), parameter :: FE_maxNnodes = 8
integer(pInt), parameter :: FE_maxNsubNodes = 56 integer(pInt), parameter :: FE_maxNsubNodes = 56
integer(pInt), parameter :: FE_maxNips = 27 integer(pInt), parameter :: FE_maxNips = 27
integer(pInt), parameter :: FE_maxNipNeighbors = 6 integer(pInt), parameter :: FE_maxNipNeighbors = 6
integer(pInt), parameter :: FE_maxmaxNnodesAtIP = 8
integer(pInt), parameter :: FE_NipFaceNodes = 4 integer(pInt), parameter :: FE_NipFaceNodes = 4
integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_Nnodes = & integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_Nnodes = &
(/8, & ! element 7 (/8, & ! element 7
@ -80,16 +82,18 @@
4, & ! element 27 4, & ! element 27
4, & ! element 157 4, & ! element 157
6, & ! element 136 6, & ! element 136
8 & ! element 21 8, & ! element 21
8 & ! element 117
/) /)
integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_NoriginalNodes = & integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_NoriginalNodes = &
(/8, & ! element 7 (/8, & ! element 7
4, & ! element 134 4, & ! element 134
4, & ! element 11 4, & ! element 11
4, & ! element 27 8, & ! element 27
4, & ! element 157 4, & ! element 157
6, & ! element 136 6, & ! element 136
20 & ! element 21 20,& ! element 21
8 & ! element 117
/) /)
integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_Nips = & integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_Nips = &
(/8, & ! element 7 (/8, & ! element 7
@ -98,7 +102,8 @@
9, & ! element 27 9, & ! element 27
4, & ! element 157 4, & ! element 157
6, & ! element 136 6, & ! element 136
27 & ! element 21 27,& ! element 21
1 & ! element 117
/) /)
integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_NipNeighbors = & integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_NipNeighbors = &
(/6, & ! element 7 (/6, & ! element 7
@ -107,7 +112,8 @@
4, & ! element 27 4, & ! element 27
6, & ! element 157 6, & ! element 157
6, & ! element 136 6, & ! element 136
6 & ! element 21 6, & ! element 21
6 & ! element 117
/) /)
integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_NsubNodes = & integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_NsubNodes = &
(/19,& ! element 7 (/19,& ! element 7
@ -116,7 +122,8 @@
12,& ! element 27 12,& ! element 27
0, & ! element 157 0, & ! element 157
15,& ! element 136 15,& ! element 136
56 & ! element 21 56,& ! element 21
0 & ! element 117
/) /)
integer(pInt), dimension(FE_maxNipNeighbors,FE_Nelemtypes), parameter :: FE_NfaceNodes = & integer(pInt), dimension(FE_maxNipNeighbors,FE_Nelemtypes), parameter :: FE_NfaceNodes = &
reshape((/& reshape((/&
@ -126,8 +133,19 @@
2,2,2,2,0,0, & ! element 27 2,2,2,2,0,0, & ! element 27
3,3,3,3,0,0, & ! element 157 3,3,3,3,0,0, & ! element 157
3,4,4,4,3,0, & ! element 136 3,4,4,4,3,0, & ! element 136
4,4,4,4,4,4 & ! element 21 4,4,4,4,4,4, & ! element 21
4,4,4,4,4,4 & ! element 117
/),(/FE_maxNipNeighbors,FE_Nelemtypes/)) /),(/FE_maxNipNeighbors,FE_Nelemtypes/))
integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_maxNnodesAtIP = &
(/1, & ! element 7
4, & ! element 134
1, & ! element 11
2, & ! element 27
1, & ! element 157
1, & ! element 136
4, & ! element 21
8 & ! element 117
/)
integer(pInt), dimension(FE_NipFaceNodes,FE_maxNipNeighbors,FE_Nelemtypes), parameter :: FE_nodeOnFace = & integer(pInt), dimension(FE_NipFaceNodes,FE_maxNipNeighbors,FE_Nelemtypes), parameter :: FE_nodeOnFace = &
reshape((/& reshape((/&
1,2,3,4 , & ! element 7 1,2,3,4 , & ! element 7
@ -171,6 +189,12 @@
3,2,6,7 , & 3,2,6,7 , &
4,3,7,8 , & 4,3,7,8 , &
4,1,5,8 , & 4,1,5,8 , &
8,7,6,5 , &
1,2,3,4 , & ! element 117
2,1,5,6 , &
3,2,6,7 , &
4,3,7,8 , &
4,1,5,8 , &
8,7,6,5 & 8,7,6,5 &
/),(/FE_NipFaceNodes,FE_maxNipNeighbors,FE_Nelemtypes/)) /),(/FE_NipFaceNodes,FE_maxNipNeighbors,FE_Nelemtypes/))
@ -277,20 +301,30 @@
integer(pInt) FE_mapElemtype integer(pInt) FE_mapElemtype
select case (what) select case (what)
case ( '7', 'C3D8') case ( '7', &
'C3D8')
FE_mapElemtype = 1 ! Three-dimensional Arbitrarily Distorted Brick FE_mapElemtype = 1 ! Three-dimensional Arbitrarily Distorted Brick
case ('134', 'C3D4') case ('134', &
'C3D4')
FE_mapElemtype = 2 ! Three-dimensional Four-node Tetrahedron FE_mapElemtype = 2 ! Three-dimensional Four-node Tetrahedron
case ( '11', 'CPE4') case ( '11', &
'CPE4')
FE_mapElemtype = 3 ! Arbitrary Quadrilateral Plane-strain FE_mapElemtype = 3 ! Arbitrary Quadrilateral Plane-strain
case ( '27', 'CPE8') case ( '27', &
'CPE8')
FE_mapElemtype = 4 ! Plane Strain, Eight-node Distorted Quadrilateral FE_mapElemtype = 4 ! Plane Strain, Eight-node Distorted Quadrilateral
case ('157') case ('157')
FE_mapElemtype = 5 ! Three-dimensional, Low-order, Tetrahedron, Herrmann Formulations FE_mapElemtype = 5 ! Three-dimensional, Low-order, Tetrahedron, Herrmann Formulations
case ('136', 'C3D6') case ('136', &
'C3D6')
FE_mapElemtype = 6 ! Three-dimensional Arbitrarily Distorted Pentahedral FE_mapElemtype = 6 ! Three-dimensional Arbitrarily Distorted Pentahedral
case ( '21', 'C3D20') case ( '21', &
FE_mapElemtype = 7 ! Three-dimensional Arbitrarily Distorted Quadratic Hexahedral 'C3D20')
FE_mapElemtype = 7 ! Three-dimensional Arbitrarily Distorted quadratic hexahedral
case ( '117', &
'123', &
'C3D8R')
FE_mapElemtype = 8 ! Three-dimensional Arbitrarily Distorted linear hexahedral with reduced integration
case default case default
FE_mapElemtype = 0 ! unknown element --> should raise an error upstream..! FE_mapElemtype = 0 ! unknown element --> should raise an error upstream..!
endselect endselect
@ -357,43 +391,50 @@
!*********************************************************** !***********************************************************
! find face-matching element of same type ! find face-matching element of same type
!*********************************************************** !***********************************************************
function mesh_faceMatch(face,elem) function mesh_faceMatch(elem,face)
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), dimension(2) :: mesh_faceMatch ! matching element's ID and corresponding face ID
integer(pInt), dimension(FE_NfaceNodes(face,mesh_element(2,elem))) :: nodeMap integer(pInt), dimension(FE_NfaceNodes(face,mesh_element(2,elem))) :: myFaceNodes ! global node ids on my face
integer(pInt) minN,NsharedElems,lonelyNode,faceNode,i,n,t integer(pInt) minN,NsharedElems, &
t, lonelyNode, &
candidateType,candidateElem, &
i,f,n
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 n = 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 myFaceNodes(n) = mesh_FEasCP('node',mesh_element(4+FE_nodeOnFace(n,face,t),elem)) ! CP id of face node
NsharedElems = mesh_sharedElem(1,nodeMap(faceNode)) ! figure # shared elements for this node NsharedElems = mesh_sharedElem(1,myFaceNodes(n)) ! figure # shared elements for this node
if (NsharedElems < minN) then if (NsharedElems < minN) then
minN = NsharedElems ! remember min # shared elems minN = NsharedElems ! remember min # shared elems
lonelyNode = faceNode ! remember most lonely node lonelyNode = n ! remember most lonely node
endif endif
enddo enddo
candidate: do i=1,minN ! iterate over lonelyNode's shared elements
mesh_faceMatch = mesh_sharedElem(1+i,nodeMap(lonelyNode)) ! present candidate elem candidate: do i = 1,minN ! iterate over lonelyNode's shared elements
if (mesh_faceMatch == elem) then ! my own element ? candidateElem = mesh_sharedElem(1+i,myFaceNodes(lonelyNode)) ! present candidate elem
mesh_faceMatch = 0_pInt ! disregard if (candidateElem == elem) cycle candidate ! my own element ?
cycle candidate candidateType = mesh_element(2,candidateElem) ! figure elemType of candidate
endif
do faceNode=1,FE_NfaceNodes(face,t) ! check remaining face nodes to match candidateFace: do f = 1,FE_maxNipNeighbors ! check each face of candidate
if (faceNode == lonelyNode) cycle ! disregard lonely node (matches anyway) if (FE_NfaceNodes(f,candidateType) /= FE_NfaceNodes(face,t)) &
n = nodeMap(faceNode) cycle candidateFace ! incompatible face
if (all(mesh_sharedElem(2:1+mesh_sharedElem(1,n),n) /= mesh_faceMatch)) then ! no ref to candidate elem? do n = 1,FE_NfaceNodes(f,candidateType) ! loop through nodes on face
mesh_faceMatch = 0_pInt ! set to "no match" (so far) if (all(myFaceNodes /= &
cycle candidate ! next candidate elem mesh_FEasCP('node', &
endif mesh_element(4+FE_nodeOnFace(n,f,candidateType),candidateElem)))) &
enddo cycle candidateFace ! other face node not one of my face nodes
exit ! surviving candidate enddo
mesh_faceMatch(1) = f
mesh_faceMatch(2) = candidateElem
exit candidate ! found my matching candidate
enddo candidateFace
enddo candidate enddo candidate
return return
@ -412,63 +453,63 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
use prec, only: pInt use prec, only: pInt
implicit none implicit none
allocate(FE_nodesAtIP(2,2,FE_maxNips,FE_Nelemtypes)) ; FE_nodesAtIP = 0_pInt allocate(FE_nodesAtIP(FE_maxmaxNnodesAtIP,FE_maxNips,FE_Nelemtypes)) ; FE_nodesAtIP = 0_pInt
allocate(FE_ipNeighbor(FE_maxNipNeighbors,FE_maxNips,FE_Nelemtypes)) ; FE_ipNeighbor = 0_pInt allocate(FE_ipNeighbor(FE_maxNipNeighbors,FE_maxNips,FE_Nelemtypes)) ; FE_ipNeighbor = 0_pInt
allocate(FE_subNodeParent(FE_maxNips,FE_maxNsubNodes,FE_Nelemtypes)) ; FE_subNodeParent = 0_pInt allocate(FE_subNodeParent(FE_maxNips,FE_maxNsubNodes,FE_Nelemtypes)) ; FE_subNodeParent = 0_pInt
allocate(FE_subNodeOnIPFace(FE_NipFaceNodes,FE_maxNipNeighbors,FE_maxNips,FE_Nelemtypes)) ; FE_subNodeOnIPFace = 0_pInt allocate(FE_subNodeOnIPFace(FE_NipFaceNodes,FE_maxNipNeighbors,FE_maxNips,FE_Nelemtypes)) ; FE_subNodeOnIPFace = 0_pInt
! fill FE_nodesAtIP with data ! fill FE_nodesAtIP with data
FE_nodesAtIP(:,:,:FE_Nips(1),1) = & ! element 7 FE_nodesAtIP(:,:FE_Nips(1),1) = & ! element 7
reshape((/& reshape((/&
1,0, 0,0, & 1, &
2,0, 0,0, & 2, &
4,0, 0,0, & 4, &
3,0, 0,0, & 3, &
5,0, 0,0, & 5, &
6,0, 0,0, & 6, &
8,0, 0,0, & 8, &
7,0, 0,0 & 7 &
/),(/2,2,FE_Nips(1)/)) /),(/FE_maxNnodesAtIP(1),FE_Nips(1)/))
FE_nodesAtIP(:,:,:FE_Nips(2),2) = & ! element 134 FE_nodesAtIP(:,:FE_Nips(2),2) = & ! element 134
reshape((/& reshape((/&
1,0, 0,0 & 1,2,3,4 &
/),(/2,2,FE_Nips(2)/)) /),(/FE_maxNnodesAtIP(2),FE_Nips(2)/))
FE_nodesAtIP(:,:,:FE_Nips(3),3) = & ! element 11 FE_nodesAtIP(:,:FE_Nips(3),3) = & ! element 11
reshape((/& reshape((/&
1,0, 0,0, & 1, &
2,0, 0,0, & 2, &
4,0, 0,0, & 4, &
3,0, 0,0 & 3 &
/),(/2,2,FE_Nips(3)/)) /),(/FE_maxNnodesAtIP(3),FE_Nips(3)/))
FE_nodesAtIP(:,:,:FE_Nips(4),4) = & ! element 27 FE_nodesAtIP(:,:FE_Nips(4),4) = & ! element 27
reshape((/& reshape((/&
1,0, 0,0, & 1,0, &
1,2, 0,0, & 1,2, &
2,0, 0,0, & 2,0, &
1,4, 0,0, & 1,4, &
1,3, 2,4, & 0,0, &
2,3, 0,0, & 2,3, &
4,0, 0,0, & 4,0, &
3,4, 0,0, & 3,4, &
3,0, 0,0 & 3,0 &
/),(/2,2,FE_Nips(4)/)) /),(/FE_maxNnodesAtIP(4),FE_Nips(4)/))
FE_nodesAtIP(:,:,:FE_Nips(5),5) = & ! element 157 FE_nodesAtIP(:,:FE_Nips(5),5) = & ! element 157
reshape((/& reshape((/&
1,0, 0,0, & 1, &
2,0, 0,0, & 2, &
3,0, 0,0, & 3, &
4,0, 0,0 & 4 &
/),(/2,2,FE_Nips(5)/)) /),(/FE_maxNnodesAtIP(5),FE_Nips(5)/))
FE_nodesAtIP(:,:,:FE_Nips(6),6) = & ! element 136 FE_nodesAtIP(:,:FE_Nips(6),6) = & ! element 136
reshape((/& reshape((/&
1,0, 0,0, & 1, &
2,0, 0,0, & 2, &
3,0, 0,0, & 3, &
4,0, 0,0, & 4, &
5,0, 0,0, & 5, &
6,0, 0,0 & 6 &
/),(/2,2,FE_Nips(6)/)) /),(/FE_maxNnodesAtIP(6),FE_Nips(6)/))
FE_nodesAtIP(:,:,:FE_Nips(7),7) = & ! element 21 FE_nodesAtIP(:,:FE_Nips(7),7) = & ! element 21
reshape((/& reshape((/&
1,0, 0,0, & 1,0, 0,0, &
1,2, 0,0, & 1,2, 0,0, &
@ -497,7 +538,7 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
8,0, 0,0, & 8,0, 0,0, &
7,8, 0,0, & 7,8, 0,0, &
7,0, 0,0 & 7,0, 0,0 &
/),(/2,2,FE_Nips(7)/)) /),(/FE_maxNnodesAtIP(7),FE_Nips(7)/))
! fill FE_ipNeighbor with data ! fill FE_ipNeighbor with data
FE_ipNeighbor(:FE_NipNeighbors(1),:FE_Nips(1),1) = & ! element 7 FE_ipNeighbor(:FE_NipNeighbors(1),:FE_Nips(1),1) = & ! element 7
@ -580,6 +621,10 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
27,25,-4,23,-6,17, & 27,25,-4,23,-6,17, &
-3,26,-4,24,-6,18 & -3,26,-4,24,-6,18 &
/),(/FE_NipNeighbors(7),FE_Nips(7)/)) /),(/FE_NipNeighbors(7),FE_Nips(7)/))
FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117
reshape((/&
-3,-5,-4,-2,-6,-1 &
/),(/FE_NipNeighbors(8),FE_Nips(8)/))
! fill FE_subNodeParent with data ! fill FE_subNodeParent with data
FE_subNodeParent(:FE_Nips(1),:FE_NsubNodes(1),1) = & ! element 7 FE_subNodeParent(:FE_Nips(1),:FE_NsubNodes(1),1) = & ! element 7
@ -604,10 +649,7 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
5, 6, 7, 8, 0, 0, 0, 0, & 5, 6, 7, 8, 0, 0, 0, 0, &
1, 2, 3, 4, 5, 6, 7, 8 & 1, 2, 3, 4, 5, 6, 7, 8 &
/),(/FE_Nips(1),FE_NsubNodes(1)/)) /),(/FE_Nips(1),FE_NsubNodes(1)/))
!FE_subNodeParent(:FE_Nips(2),:FE_NsubNodes(2),2) = & ! element 134 !FE_subNodeParent(:FE_Nips(2),:FE_NsubNodes(2),2) ! element 134 has no subnodes
! reshape((/&
! *still to be defined*
! /),(/FE_Nips(2),FE_NsubNodes(2)/))
FE_subNodeParent(:FE_Nips(3),:FE_NsubNodes(3),3) = & ! element 11 FE_subNodeParent(:FE_Nips(3),:FE_NsubNodes(3),3) = & ! element 11
reshape((/& reshape((/&
1, 2, 0, 0, & 1, 2, 0, 0, &
@ -712,6 +754,7 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
7, 7, 7, 7, 7, 7, 7, 7, 3, 3, 3, 3, 6, 6, 6, 6, 8, 8, 8, 8, 2, 2, 4, 4, 5, 5, 1, & 7, 7, 7, 7, 7, 7, 7, 7, 3, 3, 3, 3, 6, 6, 6, 6, 8, 8, 8, 8, 2, 2, 4, 4, 5, 5, 1, &
8, 8, 8, 8, 8, 8, 8, 8, 4, 4, 4, 4, 5, 5, 5, 5, 7, 7, 7, 7, 1, 1, 3, 3, 6, 6, 2 & 8, 8, 8, 8, 8, 8, 8, 8, 4, 4, 4, 4, 5, 5, 5, 5, 7, 7, 7, 7, 1, 1, 3, 3, 6, 6, 2 &
/),(/FE_Nips(7),FE_NsubNodes(7)/)) /),(/FE_Nips(7),FE_NsubNodes(7)/))
!FE_subNodeParent(:FE_Nips(8),:FE_NsubNodes(8),8) ! element 117 has no subnodes
! fill FE_subNodeOnIPFace with data ! fill FE_subNodeOnIPFace with data
FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(1),:FE_Nips(1),1) = & ! element 7 FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(1),:FE_Nips(1),1) = & ! element 7
@ -1038,6 +1081,15 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
55,28, 7,29, & 55,28, 7,29, &
63,49,23,48 & 63,49,23,48 &
/),(/FE_NipFaceNodes,FE_NipNeighbors(7),FE_Nips(7)/)) /),(/FE_NipFaceNodes,FE_NipNeighbors(7),FE_Nips(7)/))
FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117
reshape((/&
2, 3, 7, 6, & ! 1
1, 5, 8, 4, &
3, 4, 8, 7, &
1, 2, 6, 5, &
5, 6, 7, 8, &
1, 4, 3, 2 &
/),(/FE_NipFaceNodes,FE_NipNeighbors(8),FE_Nips(8)/))
return return
@ -1111,18 +1163,18 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
pos = IO_stringPos(line,maxNchunks) pos = IO_stringPos(line,maxNchunks)
if ( IO_lc(IO_StringValue(line,pos,1)) == 'resolution') then if ( IO_lc(IO_StringValue(line,pos,1)) == 'resolution') then
do i = 2,6,2 do i = 2,6,2
select case (IO_lc(IO_stringValue(line,pos,i))) select case (IO_lc(IO_stringValue(line,pos,i)))
case('a') case('a')
a = IO_intValue(line,pos,i+1) a = IO_intValue(line,pos,i+1)
case('b') case('b')
b = IO_intValue(line,pos,i+1) b = IO_intValue(line,pos,i+1)
case('c') case('c')
c = IO_intValue(line,pos,i+1) c = IO_intValue(line,pos,i+1)
end select end select
enddo enddo
mesh_Nelems = 2**(a+b+c) mesh_Nelems = 2**(a+b+c)
mesh_Nnodes = (1+2**a)*(1+2**b)*(1+2**c) mesh_Nnodes = (1+2**a)*(1+2**b)*(1+2**c)
exit exit
endif endif
enddo enddo
@ -2055,31 +2107,31 @@ subroutine mesh_marc_count_cpSizes (unit)
select case ( IO_lc(IO_StringValue(line,pos,1)) ) select case ( IO_lc(IO_StringValue(line,pos,1)) )
case ('resolution') case ('resolution')
gotResolution = .true. gotResolution = .true.
do i = 2,6,2 do i = 2,6,2
tag = IO_lc(IO_stringValue(line,pos,i)) tag = IO_lc(IO_stringValue(line,pos,i))
select case (tag) select case (tag)
case('a') case('a')
a = 1+2**IO_intValue(line,pos,i+1) a = 1+2**IO_intValue(line,pos,i+1)
case('b') case('b')
b = 1+2**IO_intValue(line,pos,i+1) b = 1+2**IO_intValue(line,pos,i+1)
case('c') case('c')
c = 1+2**IO_intValue(line,pos,i+1) c = 1+2**IO_intValue(line,pos,i+1)
end select end select
enddo enddo
case ('dimension') case ('dimension')
gotDimension = .true. gotDimension = .true.
do i = 2,6,2 do i = 2,6,2
tag = IO_lc(IO_stringValue(line,pos,i)) tag = IO_lc(IO_stringValue(line,pos,i))
select case (tag) select case (tag)
case('x') case('x')
x = IO_floatValue(line,pos,i+1) x = IO_floatValue(line,pos,i+1)
case('y') case('y')
y = IO_floatValue(line,pos,i+1) y = IO_floatValue(line,pos,i+1)
case('z') case('z')
z = IO_floatValue(line,pos,i+1) z = IO_floatValue(line,pos,i+1)
end select end select
enddo enddo
end select end select
if (gotDimension .and. gotResolution) exit if (gotDimension .and. gotResolution) exit
enddo enddo
@ -2240,21 +2292,21 @@ subroutine mesh_marc_count_cpSizes (unit)
case ('homogenization') case ('homogenization')
gotHomogenization = .true. gotHomogenization = .true.
homog = IO_intValue(line,pos,2) homog = IO_intValue(line,pos,2)
case ('resolution') case ('resolution')
gotResolution = .true. gotResolution = .true.
do i = 2,6,2 do i = 2,6,2
select case (IO_lc(IO_stringValue(line,pos,i))) select case (IO_lc(IO_stringValue(line,pos,i)))
case('a') case('a')
a = 2**IO_intValue(line,pos,i+1) a = 2**IO_intValue(line,pos,i+1)
case('b') case('b')
b = 2**IO_intValue(line,pos,i+1) b = 2**IO_intValue(line,pos,i+1)
case('c') case('c')
c = 2**IO_intValue(line,pos,i+1) c = 2**IO_intValue(line,pos,i+1)
end select end select
enddo enddo
end select end select
if (gotDimension .and. gotHomogenization .and. gotResolution) exit if (gotDimension .and. gotHomogenization .and. gotResolution) exit
enddo enddo
@ -2267,18 +2319,18 @@ subroutine mesh_marc_count_cpSizes (unit)
pos = IO_stringPos(line,1) pos = IO_stringPos(line,1)
e = e+1 ! valid element entry e = e+1 ! valid element entry
mesh_element ( 1,e) = e ! FE id mesh_element ( 1,e) = e ! FE id
mesh_element ( 2,e) = FE_mapElemtype('C3D8R') ! elem type mesh_element ( 2,e) = FE_mapElemtype('C3D8R') ! elem type
mesh_element ( 3,e) = homog ! homogenization mesh_element ( 3,e) = homog ! homogenization
mesh_element ( 4,e) = IO_IntValue(line,pos,1) ! microstructure mesh_element ( 4,e) = IO_IntValue(line,pos,1) ! microstructure
mesh_element ( 5,e) = (e-1) + (e-1)/a + (e-1)/a/b*(a+1) ! base node mesh_element ( 5,e) = (e-1) + (e-1)/a + (e-1)/a/b*(a+1) ! base node
mesh_element ( 6,e) = mesh_element ( 5,e) + 1 mesh_element ( 6,e) = mesh_element ( 5,e) + 1
mesh_element ( 7,e) = mesh_element ( 5,e) + (a+1) + 1 mesh_element ( 7,e) = mesh_element ( 5,e) + (a+1) + 1
mesh_element ( 8,e) = mesh_element ( 5,e) + (a+1) mesh_element ( 8,e) = mesh_element ( 5,e) + (a+1)
mesh_element ( 9,e) = mesh_element ( 5,e) + (a+1)*(b+1) ! second floor base node mesh_element ( 9,e) = mesh_element ( 5,e) + (a+1)*(b+1) ! second floor base node
mesh_element (10,e) = mesh_element ( 9,e) + 1 mesh_element (10,e) = mesh_element ( 9,e) + 1
mesh_element (11,e) = mesh_element ( 9,e) + (a+1) + 1 mesh_element (11,e) = mesh_element ( 9,e) + (a+1) + 1
mesh_element (12,e) = mesh_element ( 9,e) + (a+1) mesh_element (12,e) = mesh_element ( 9,e) + (a+1)
enddo enddo
110 return 110 return
@ -2545,9 +2597,13 @@ subroutine mesh_marc_count_cpSizes (unit)
use prec, only: pInt use prec, only: pInt
implicit none implicit none
integer(pInt) e,t,i,j,k,n integer(pInt) e,t,i,j,k,l,m,n,a,anchor, neighborType
integer(pInt) neighbor,neighboringElem,neighboringIP,matchingElem integer(pInt) neighbor,neighboringElem,neighboringIP
integer(pInt), dimension(2) :: linkedNode = 0_pInt integer(pInt), dimension(2) :: matchingElem ! face and id of matching elem
integer(pInt), dimension(FE_maxmaxNnodesAtIP) :: linkedNodes, &
matchingNodes
linkedNodes = 0_pInt
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
@ -2562,42 +2618,39 @@ subroutine mesh_marc_count_cpSizes (unit)
else ! neighboring element's IP else ! neighboring element's IP
neighboringElem = 0_pInt neighboringElem = 0_pInt
neighboringIP = 0_pInt neighboringIP = 0_pInt
matchingElem = mesh_faceMatch(-neighbor,e) ! get CP elem id of face match matchingElem = mesh_faceMatch(e,-neighbor) ! get face and CP elem id of face match
if (matchingElem > 0 .and. mesh_element(2,matchingElem) == t) then ! found match of same type? if (matchingElem(2) > 0_pInt) then ! found match?
if (FE_nodesAtIP(2,1,i,t) == 0) then ! single linked node neighborType = mesh_element(2,matchingElem(2))
matchNode1: do j = 1,FE_Nnodes(t) ! check against all neighbor's nodes l = 0_pInt
if (mesh_element(4+FE_nodesAtIP(1,1,i,t),e)==mesh_element(4+j,matchingElem)) then do a = 1,FE_maxNnodesAtIP(t) ! figure my anchor nodes on connecting face
linkedNode(1) = j ! which neighboring node matches my first nodeAtIP (indexed globally) anchor = FE_nodesAtIP(a,i,t)
linkedNode(2) = 0_pInt if (any(FE_nodeOnFace(:,-neighbor,t) == anchor)) then ! ip anchor sits on face?
exit matchNode1 l = l+1_pInt
linkedNodes(l) = mesh_element(4+anchor,e) ! FE id of anchor node
endif
enddo
neighborIP: do j = 1,FE_Nips(neighborType) ! loop over neighboring ips
m = 0_pInt
matchingNodes = 0_pInt
do a = 1,FE_maxNnodesAtIP(neighborType) ! check each anchor node of that ip
anchor = FE_nodesAtIP(a,j,neighborType)
if ( anchor /= 0_pInt .and. & ! valid anchor node
any(FE_nodeOnFace(:,matchingElem(1),neighborType) == anchor) ) then
m = m+1_pInt
matchingNodes(m) = mesh_element(4+anchor,matchingElem(2)) ! FE id of neighbor's anchor node
endif endif
enddo matchNode1 enddo
matchFace1: do j = 1,FE_Nips(t) if (m /= l) cycle neighborIP ! this ip has wrong count of anchors on face
if ( (linkedNode(1) == FE_nodesAtIP(1,1,j,t)) .and. & do a = 1,l
(FE_nodesAtIP(2,1,j,t) == 0) ) then if (all(matchingNodes /= linkedNodes(a))) cycle neighborIP
neighboringElem = matchingElem enddo
neighboringIP = j neighboringElem = matchingElem(2) ! survival of the fittest
exit matchFace1 neighboringIP = j
endif exit neighborIP
enddo matchFace1 enddo neighborIP
else ! double linked node endif ! end of valid external matching
matchNode2: do j = 1,FE_Nnodes(t) ! check against all neighbor's nodes endif ! end of internal/external matching
if (mesh_element(4+FE_nodesAtIP(1,1,i,t),e)==mesh_element(4+j,matchingElem)) linkedNode(1) = j ! which neighboring node matches my first nodeAtIP (indexed globally)
if (mesh_element(4+FE_nodesAtIP(2,1,i,t),e)==mesh_element(4+j,matchingElem)) linkedNode(2) = j ! which neighboring node matches my second nodeAtIP (indexed globally)
enddo matchNode2
matchFace2: do j = 1,FE_Nips(t)
if ((linkedNode(1) == FE_nodesAtIP(1,1,j,t) .and. linkedNode(2) == FE_nodesAtIP(2,1,j,t)) .or. &
(linkedNode(1) == FE_nodesAtIP(2,1,j,t) .and. linkedNode(2) == FE_nodesAtIP(1,1,j,t)) .or. &
(linkedNode(1) == FE_nodesAtIP(1,2,j,t) .and. linkedNode(2) == FE_nodesAtIP(2,2,j,t)) .or. &
(linkedNode(1) == FE_nodesAtIP(2,2,j,t) .and. linkedNode(2) == FE_nodesAtIP(1,2,j,t))) then
neighboringElem = matchingElem
neighboringIP = j
exit matchFace2
endif
enddo matchFace2
endif
endif
endif
mesh_ipNeighborhood(1,n,i,e) = neighboringElem mesh_ipNeighborhood(1,n,i,e) = neighboringElem
mesh_ipNeighborhood(2,n,i,e) = neighboringIP mesh_ipNeighborhood(2,n,i,e) = neighboringIP
enddo enddo
@ -2737,17 +2790,17 @@ subroutine mesh_marc_count_cpSizes (unit)
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 interfaces of IP do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP
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)
forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles) ! start at each interface node and build valid triangles to cover interface forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles) ! start at each interface node and build valid triangles to cover interface
normal(:,j,n) = math_vectorproduct(nPos(:,1+mod(n+j-1,FE_NipFaceNodes)) - nPos(:,n), & ! calc their normal vectors normal(:,j,n) = math_vectorproduct(nPos(:,1+mod(n+j-1,FE_NipFaceNodes)) - nPos(:,n), & ! calc their normal vectors
nPos(:,1+mod(n+j-0,FE_NipFaceNodes)) - nPos(:,n)) nPos(:,1+mod(n+j-0,FE_NipFaceNodes)) - nPos(:,n))
area(j,n) = dsqrt(sum(normal(:,j,n)*normal(:,j,n))) ! and area area(j,n) = dsqrt(sum(normal(:,j,n)*normal(:,j,n))) ! and area
end forall end forall
forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles, area(j,n) > 0.0_pReal) & forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles, area(j,n) > 0.0_pReal) &
normal(:,j,n) = normal(:,j,n) / area(j,n) ! make unit normal normal(:,j,n) = normal(:,j,n) / area(j,n) ! make unit normal
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) / (FE_NipFaceNodes*2.0_pReal) ! area of parallelograms instead of triangles
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(sum(normal,3),2) / count(area > 0.0_pReal) ! average of all valid normals
enddo enddo
enddo enddo
enddo enddo
@ -2783,7 +2836,7 @@ subroutine mesh_marc_count_cpSizes (unit)
mesh_HomogMicro(mesh_element(3,i),mesh_element(4,i)) + 1 ! count combinations of homogenization and microstructure mesh_HomogMicro(mesh_element(3,i),mesh_element(4,i)) + 1 ! count combinations of homogenization and microstructure
enddo enddo
!$OMP CRITICAL (write2out) $OMP CRITICAL (write2out)
! write(6,*) ! write(6,*)
! write(6,*) 'Input Parser: IP COORDINATES' ! write(6,*) 'Input Parser: IP COORDINATES'