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

Simplified 'mesh_get_nodeElemDimensions', mesh_build_elements' and 'mesh_build_sharedElems'. 

Corrected entries in 'FE_subNodeOnIPFace'
This commit is contained in:
Christoph Kords 2009-04-03 07:03:25 +00:00
parent b2e8510b6d
commit 0e84b07dc7
1 changed files with 218 additions and 244 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

462
trunk/mesh.f90
View File

@ -26,7 +26,8 @@
! FE-solver specific part (different for MARC/ABAQUS)..! ! FE-solver specific part (different for MARC/ABAQUS)..!
! Hence, I suggest to prefix with "FE_" ! Hence, I suggest to prefix with "FE_"
! !
! _Nnodes : # nodes in a specific type of element ! _Nnodes : # nodes in a specific type of element (how we use it)
! _NoriginalNodes : # nodes in a specific type of element (how it is originally defined by marc)
! _Nips : # IPs in a specific type of element ! _Nips : # IPs in a specific type of element
! _NipNeighbors : # IP neighbors in a specific type of element ! _NipNeighbors : # IP neighbors in a specific type of element
! _ipNeighbor : +x,-x,+y,-y,+z,-z list of intra-element IPs and ! _ipNeighbor : +x,-x,+y,-y,+z,-z list of intra-element IPs and
@ -58,12 +59,21 @@
integer(pInt) :: hypoelasticTableStyle = 0 integer(pInt) :: hypoelasticTableStyle = 0
integer(pInt) :: initialcondTableStyle = 0 integer(pInt) :: initialcondTableStyle = 0
integer(pInt), parameter :: FE_Nelemtypes = 7 integer(pInt), parameter :: FE_Nelemtypes = 7
integer(pInt), parameter :: FE_maxNnodes = 20 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_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
4, & ! element 134
4, & ! element 11
4, & ! element 27
4, & ! element 157
6, & ! element 136
8 & ! element 21
/)
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
@ -929,8 +939,8 @@
3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 4, 4, 4, 4, 7, 7, 7, 7, 1, 1, 6, 6, 8, 8, 5, & 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 4, 4, 4, 4, 7, 7, 7, 7, 1, 1, 6, 6, 8, 8, 5, &
4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 1, 3, 3, 3, 3, 8, 8, 8, 8, 2, 2, 5, 5, 7, 7, 6, & 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 1, 3, 3, 3, 3, 8, 8, 8, 8, 2, 2, 5, 5, 7, 7, 6, &
5, 5, 5, 5, 5, 5, 5, 5, 1, 1, 1, 1, 6, 6, 6, 6, 8, 8, 8, 8, 2, 2, 4, 4, 7, 7, 3, & 5, 5, 5, 5, 5, 5, 5, 5, 1, 1, 1, 1, 6, 6, 6, 6, 8, 8, 8, 8, 2, 2, 4, 4, 7, 7, 3, &
6, 6, 6, 6, 6, 6, 6, 6, 2, 2, 2, 2, 5, 5, 5, 5, 8, 8, 8, 8, 1, 1, 3, 3, 8, 8, 4, & 6, 6, 6, 6, 6, 6, 6, 6, 2, 2, 2, 2, 5, 5, 5, 5, 7, 7, 7, 7, 1, 1, 3, 3, 8, 8, 4, &
7, 7, 7, 7, 7, 7, 7, 7, 3, 3, 3, 3, 6, 6, 6, 6, 7, 7, 7, 7, 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_maxNips,FE_maxNsubNodes,FE_Nelemtypes/)) /),(/FE_maxNips,FE_maxNsubNodes,FE_Nelemtypes/))
integer(pInt), dimension(FE_NipFaceNodes,FE_maxNipNeighbors,FE_maxNips,FE_Nelemtypes), parameter :: FE_subNodeOnIPFace = & integer(pInt), dimension(FE_NipFaceNodes,FE_maxNipNeighbors,FE_maxNips,FE_Nelemtypes), parameter :: FE_subNodeOnIPFace = &
@ -1907,168 +1917,168 @@
0, 0, 0, 0, & 0, 0, 0, 0, &
0, 0, 0, 0, & 0, 0, 0, 0, &
0, 0, 0, 0, & 0, 0, 0, 0, &
9,33,57,37, & ! element 21 9,33,57,37, & ! element 21
1,17,44,16, & 1,17,44,16, &
16,44,57,33, & 33,16,44,57, &
1, 9,37,17, & 1, 9,37,17, &
17,37,44,57, & 17,37,57,44, &
1,16,33, 9, & 1,16,33, 9, &
10,34,58,38, & ! 10,34,58,38, & ! 2
9,37,57,33, & 9,37,57,33, &
33,57,58,34, & 34,33,57,58, &
9,10,38,37, & 9,10,38,37, &
37,38,57,58, & 37,38,58,57, &
9,33,34,10, & 9,33,34,10, &
2,11,39,18, & ! 2,11,39,18, & ! 3
10,38,58,34, & 10,38,58,34, &
34,58,39,11, & 11,34,58,39, &
10, 2,18,38, & 10, 2,18,38, &
38,18,58,39, & 38,18,39,58, &
10,34,11, 2, & 10,34,11, 2, &
16,15,43,44, & ! 33,36,60,57, & ! 4
2,18,39,11, &
11,39,43,15, &
2,16,44,18, &
18,44,39,43, &
2,11,15,16, &
33,36,60,57, & !
16,44,43,15, & 16,44,43,15, &
15,43,60,36, & 36,15,43,60, &
16,33,57,44, & 16,33,57,44, &
44,57,43,60, & 44,57,60,43, &
16,15,36,33, & 16,15,36,33, &
34,35,59,58, & ! 34,35,59,58, & ! 5
33,57,60,36, & 33,57,60,36, &
36,60,59,35, & 35,36,60,59, &
33,34,58,57, & 33,34,58,57, &
57,58,60,59, & 57,58,59,60, &
33,36,35,34, & 33,36,35,34, &
11,12,40,39, & ! 11,12,40,39, & ! 6
34,58,59,35, & 34,58,59,35, &
35,59,40,12, & 12,35,59,40, &
34,11,39,58, & 34,11,39,58, &
58,39,59,40, & 58,39,40,59, &
34,35,12,11, & 34,35,12,11, &
15, 4,20,43, & ! 36,14,42,60, & ! 7
11,39,40,12, &
12,40,20, 4, &
11,15,43,39, &
39,43,40,20, &
11,12, 4,15, &
36,14,42,60, & !
15,43,20, 4, & 15,43,20, 4, &
4,20,42,14, & 14, 4,20,42, &
15,36,60,43, & 15,36,60,43, &
43,60,20,42, & 43,60,42,20, &
15, 4,14,36, & 15, 4,14,36, &
37,57,61,45, & ! 35,13,41,59, & ! 8
36,60,42,14, &
13,14,42,41, &
36,35,59,60, &
60,59,41,42, &
36,14,13,35, &
12, 3,19,40, & ! 9
35,59,41,13, &
3,13,41,19, &
35,12,40,59, &
59,40,19,41, &
35,13, 3,12, &
37,57,61,45, & ! 10
17,21,52,44, & 17,21,52,44, &
44,52,61,57, & 57,44,52,61, &
17,37,45,21, & 17,37,45,21, &
21,45,52,61, & 21,45,61,52, &
17,44,57,37, & 17,44,57,37, &
38,58,62,46, & ! 38,58,62,46, & ! 11
37,45,61,57, & 37,45,61,57, &
57,61,62,58, & 58,57,61,62, &
37,38,46,45, & 37,38,46,45, &
45,46,61,62, & 45,46,62,61, &
37,57,58,38, & 37,57,58,38, &
18,39,47,22, & ! 18,39,47,22, & ! 12
38,46,62,58, & 38,46,62,58, &
58,62,47,39, & 39,58,62,47, &
38,18,22,46, & 38,18,22,46, &
46,22,62,47, & 46,22,47,62, &
38,58,39,18, & 38,58,39,18, &
44,43,51,52, & ! 57,60,64,61, & ! 13
18,22,47,39, &
39,47,51,43, &
18,44,52,22, &
22,52,47,51, &
18,39,43,44, &
57,60,64,61, & !
44,52,51,43, & 44,52,51,43, &
43,51,64,60, & 60,43,51,64, &
44,57,61,52, & 44,57,61,52, &
52,61,51,64, & 52,61,64,51, &
44,43,60,57, & 44,43,60,57, &
58,59,63,62, & ! 58,59,63,62, & ! 14
57,61,64,60, & 57,61,64,60, &
60,64,63,59, & 59,60,64,63, &
57,58,62,61, & 57,58,62,61, &
61,62,64,63, & 61,62,63,64, &
57,60,59,58, & 57,60,59,58, &
39,40,48,47, & ! 39,40,48,47, & ! 15
58,62,63,59, & 58,62,63,59, &
59,63,48,40, & 40,59,63,48, &
58,39,47,62, & 58,39,47,62, &
62,47,63,48, & 62,47,48,63, &
58,59,40,39, & 58,59,40,39, &
43,20,24,51, & ! 60,42,50,64, & ! 16
39,47,48,40, &
40,48,24,20, &
39,43,51,47, &
47,51,48,24, &
39,40,20,43, &
60,42,50,64, & !
43,51,24,20, & 43,51,24,20, &
20,24,50,42, & 42,20,24,50, &
43,60,64,51, & 43,60,64,51, &
51,64,24,50, & 51,64,50,24, &
43,20,42,60, & 43,20,42,60, &
45,61,53,25, & ! 59,41,49,63, & ! 17
60,64,50,42, &
41,42,50,49, &
60,59,63,64, &
64,63,49,50, &
60,42,41,59, &
40,19,23,48, & ! 18
59,63,49,41, &
19,41,49,23, &
59,40,48,63, &
63,48,23,49, &
59,41,19,40, &
45,61,53,25, & ! 19
21, 5,32,52, & 21, 5,32,52, &
52,32,53,61, & 61,52,32,53, &
21,45,25, 5, & 21,45,25, 5, &
5,25,32,53, & 5,25,53,32, &
21,52,61,45, & 21,52,61,45, &
46,62,54,26, & ! 46,62,54,26, & ! 20
45,25,53,61, & 45,25,53,61, &
61,53,54,62, & 62,61,53,54, &
45,46,26,25, & 45,46,26,25, &
25,26,53,54, & 25,26,54,53, &
45,61,62,46, & 45,61,62,46, &
22,47,27, 6, & ! 22,47,27, 6, & ! 21
46,26,54,62, & 46,26,54,62, &
62,54,27,47, & 47,62,54,27, &
46,22, 6,26, & 46,22, 6,26, &
26, 6,54,27, & 26, 6,27,54, &
46,62,47,22, & 46,62,47,22, &
52,51,31,32, & ! 61,64,56,53, & ! 22
22, 6,27,47, &
47,27,31,51, &
22,52,32, 6, &
6,32,27,31, &
22,47,51,52, &
61,64,56,53, & !
52,32,31,51, & 52,32,31,51, &
51,31,56,64, & 64,51,31,56, &
52,61,53,32, & 52,61,53,32, &
32,53,31,56, & 32,53,56,31, &
52,51,64,61, & 52,51,64,61, &
62,63,55,54, & ! 62,63,55,54, & ! 23
61,53,56,64, & 61,53,56,64, &
64,56,55,63, & 63,64,56,55, &
61,62,54,53, & 61,62,54,53, &
53,54,56,55, & 53,54,55,56, &
61,64,63,62, & 61,64,63,62, &
47,48,28,27, & ! 47,48,28,27, & ! 24
62,54,55,63, & 62,54,55,63, &
63,55,28,48, & 48,63,55,28, &
62,47,27,54, & 62,47,27,54, &
54,27,55,28, & 54,27,28,55, &
62,63,48,47, & 62,63,48,47, &
51,24, 8,31, & ! 64,50,30,56, & ! 25
47,27,28,48, &
48,28, 8,24, &
47,51,31,27, &
27,31,28, 8, &
47,48,24,51, &
64,50,30,56, & !
51,31, 8,24, & 51,31, 8,24, &
24, 8,30,50, & 50,24, 8,30, &
51,64,56,31, & 51,64,56,31, &
31,56, 8,30, & 31,56,30, 8, &
51,24,50,64 & 51,24,50,64, &
63,49,29,55, & ! 26
64,56,30,50, &
49,50,30,29, &
64,63,55,56, &
56,55,29,30, &
64,50,49,63, &
48,23, 7,28, & ! 27
63,55,29,49, &
23,49,29, 7, &
63,48,28,55, &
55,28, 7,29, &
63,49,23,48 &
/),(/FE_NipFaceNodes,FE_maxNipNeighbors,FE_maxNips,FE_Nelemtypes/)) /),(/FE_NipFaceNodes,FE_maxNipNeighbors,FE_maxNips,FE_Nelemtypes/))
CONTAINS CONTAINS
@ -2326,16 +2336,17 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
! assign globals: ! assign globals:
! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsharedElems ! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsharedElems
!******************************************************************** !********************************************************************
SUBROUTINE mesh_get_nodeElemDimensions (unit) SUBROUTINE mesh_get_nodeElemDimensions (unit)
use prec, only: pInt use prec, only: pInt
use IO use IO
implicit none implicit none
integer(pInt), parameter :: maxNchunks = 66
integer(pInt), dimension (mesh_Nnodes) :: node_count integer(pInt), dimension (mesh_Nnodes) :: node_count
integer(pInt), dimension (:), allocatable :: node_seen integer(pInt), dimension (:), allocatable :: node_seen
integer(pInt) unit,i,j,n,t,e,NnodesStillToBeRead integer(pInt) unit,i,j,n,t,e
integer(pInt), dimension (133) :: pos integer(pInt), dimension (1+2*maxNchunks) :: pos
character*300 line character*300 line
610 FORMAT(A300) 610 FORMAT(A300)
@ -2348,50 +2359,35 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
read (unit,610,END=630) line read (unit,610,END=630) line
pos = IO_stringPos(line,1) pos = IO_stringPos(line,1)
if( IO_lc(IO_stringValue(line,pos,1)) == 'connectivity' ) then if( IO_lc(IO_stringValue(line,pos,1)) == 'connectivity' ) then
read (unit,610,END=630) line ! Garbage line read (unit,610,END=630) line ! Garbage line
readThisElement: do i=1,mesh_Nelems ! read all elements do i=1,mesh_Nelems ! read all elements
NnodesStillToBeRead = -1_pInt ! number of nodes per element still to be read read (unit,610,END=630) line
do while (NnodesStillToBeRead /= 0) pos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max (plus ID, type)
read (unit,610,END=630) line e = mesh_FEasCP('elem',IO_intValue(line,pos,1))
pos = IO_stringPos(line,2+FE_maxNnodes) ! limit to maximal number of nodes (plus ID, type) if (e /= 0) then
if (NnodesStillToBeRead < 0) then ! first line of element with FE ID and element type t = FE_mapElemtype(IO_StringValue(line,pos,2))
e = mesh_FEasCP('elem',IO_intValue(line,pos,1)) mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(t))
if (e < 1) exit readThisElement ! disregard non-CP elems mesh_maxNips = max(mesh_maxNips,FE_Nips(t))
t = FE_mapElemtype(IO_StringValue(line,pos,2)) ! element type mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(t))
NnodesStillToBeRead = FE_Nnodes(t) mesh_maxNsubNodes = max(mesh_maxNsubNodes,FE_NsubNodes(t))
mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(t)) node_seen = 0_pInt
mesh_maxNips = max(mesh_maxNips,FE_Nips(t)) do j=1,FE_Nnodes(t)
mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(t)) n = mesh_FEasCP('node',IO_IntValue (line,pos,j+2))
mesh_maxNsubNodes = max(mesh_maxNsubNodes,FE_NsubNodes(t)) if (all(node_seen /= n)) node_count(n) = node_count(n)+1
node_seen = 0_pInt node_seen(j) = n
do j = 1,(pos(1)-2) end do
n = mesh_FEasCP('node',IO_IntValue (line,pos,j+2)) call IO_skipChunks(unit,FE_NoriginalNodes(t)-(pos(1)-2)) ! read on if FE_Nnodes exceeds node count present on current line
if (all(node_seen /= n)) then end if
node_count(n) = node_count(n)+1 end do
end if
node_seen(j) = n
NnodesStillToBeRead = NnodesStillToBeRead - 1
end do
else ! all lines after the first just contain nodes
do j = 1,pos(1)
n = mesh_FEasCP('node',IO_IntValue (line,pos,j))
if (all(node_seen /= n)) then
node_count(n) = node_count(n)+1
end if
node_seen(FE_Nnodes(t)-NnodesStillToBeRead+1) = n
NnodesStillToBeRead = NnodesStillToBeRead - 1
end do
end if
end do
end do readThisElement
exit exit
end if end if
end do end do
630 mesh_maxNsharedElems = maxval(node_count) 630 mesh_maxNsharedElems = maxval(node_count)
return return
END SUBROUTINE END SUBROUTINE
!******************************************************************** !********************************************************************
! Build element set mapping ! Build element set mapping
! !
@ -2442,9 +2438,10 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
use IO use IO
implicit none implicit none
integer(pInt), parameter :: maxNchunks = 1
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 (1+2*maxNchunks) :: pos
character*300 line character*300 line
610 FORMAT(A300) 610 FORMAT(A300)
@ -2455,7 +2452,7 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
rewind(unit) rewind(unit)
do do
read (unit,610,END=650) line read (unit,610,END=650) line
pos = IO_stringPos(line,1) pos = IO_stringPos(line,maxNchunks)
if( IO_lc(IO_stringValue(line,pos,1)) == 'coordinates' ) then if( IO_lc(IO_stringValue(line,pos,1)) == 'coordinates' ) then
read (unit,610,END=650) line ! skip crap line read (unit,610,END=650) line ! skip crap line
do i=1,mesh_Nnodes do i=1,mesh_Nnodes
@ -2577,8 +2574,9 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
use IO use IO
implicit none implicit none
integer unit,e,i,j,sv,val,NnodesStillToBeRead integer unit,i,j,sv,val,e
integer, parameter :: maxNchunks = 2+FE_maxNnodes
integer(pInt), parameter :: maxNchunks = 66
integer(pInt), dimension(1+2*maxNchunks) :: pos integer(pInt), dimension(1+2*maxNchunks) :: pos
integer(pInt), dimension(1+mesh_NcpElems) :: contInts integer(pInt), dimension(1+mesh_NcpElems) :: contInts
character*300 line character*300 line
@ -2587,76 +2585,66 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
610 FORMAT(A300) 610 FORMAT(A300)
rewind(unit) rewind(unit)
do do
read (unit,610,END=680) line read (unit,610,END=620) line
pos = IO_stringPos(line,2) pos = IO_stringPos(line,2)
if( IO_lc(IO_stringValue(line,pos,1)) == 'connectivity' ) then if( IO_lc(IO_stringValue(line,pos,1)) == 'connectivity' ) then
read (unit,610,END=680) line ! Garbage line read (unit,610,END=620) line ! Garbage line
readThisElement: do i=1,mesh_Nelems do i=1,mesh_Nelems
NnodesStillToBeRead = -1_pInt ! number of nodes per element still to be read read (unit,610,END=620) line
do while (NnodesStillToBeRead /= 0) pos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max (plus ID, type)
read (unit,610,END=680) line e = mesh_FEasCP('elem',IO_intValue(line,pos,1))
pos = IO_stringPos(line,2+FE_maxNnodes) ! limit to maximal number of nodes (plus ID, type) if (e /= 0) then ! disregard non CP elems
if (NnodesStillToBeRead < 0) then ! first line of element with FE ID and element type mesh_element (1,e) = IO_IntValue (line,pos,1) ! FE id
e = mesh_FEasCP('elem',IO_intValue(line,pos,1)) mesh_element (2,e) = FE_mapElemtype(IO_StringValue (line,pos,2)) ! elem type
if (e < 1) exit readThisElement ! disregard non-CP elems forall (j=1:FE_Nnodes(mesh_element(2,e))) &
mesh_element(1,e) = IO_IntValue(line,pos,1) ! FE id mesh_element(j+4,e) = IO_IntValue (line,pos,j+2) ! copy FE ids of nodes
mesh_element(2,e) = FE_mapElemtype(IO_StringValue(line,pos,2)) ! elem type call IO_skipChunks(unit,FE_NoriginalNodes(mesh_element(2,e))-(pos(1)-2)) ! read on if FE_Nnodes exceeds node count present on current line
NnodesStillToBeRead = FE_Nnodes(mesh_element(2,e)) end if
do j = 1,(pos(1)-2) end do
mesh_element(j+4,e) = IO_IntValue(line,pos,j+2) ! copy FE ids of nodes
NnodesStillToBeRead = NnodesStillToBeRead - 1
end do
else ! all lines after the first just contain nodes
do j = 1,pos(1)
mesh_element(FE_Nnodes(mesh_element(2,e))-NnodesStillToBeRead+5,e) = IO_IntValue(line,pos,j) ! copy FE ids of nodes
NnodesStillToBeRead = NnodesStillToBeRead - 1
end do
end if
end do
end do readThisElement
exit exit
end if endif
end do enddo
rewind(unit) ! just in case "initial state" apears before "connectivity" rewind(unit) ! just in case "initial state" apears before "connectivity"
read (unit,610,END=680) line read (unit,610,END=620) line
do do
pos = IO_stringPos(line,2) pos = IO_stringPos(line,2)
if( (IO_lc(IO_stringValue(line,pos,1)) == 'initial').and. & if( (IO_lc(IO_stringValue(line,pos,1)) == 'initial').and. &
(IO_lc(IO_stringValue(line,pos,2)) == 'state') ) then (IO_lc(IO_stringValue(line,pos,2)) == 'state') ) then
if (initialcondTableStyle == 2) read (unit,610,END=680) line ! read extra line for new style if (initialcondTableStyle == 2) read (unit,610,END=620) line ! read extra line for new style
read (unit,610,END=680) line ! read line with index of state var read (unit,610,END=620) line ! read line with index of state var
pos = IO_stringPos(line,1) pos = IO_stringPos(line,1)
sv = IO_IntValue (line,pos,1) ! figure state variable index sv = IO_IntValue (line,pos,1) ! figure state variable index
if( (sv == 2).or.(sv == 3) ) then ! only state vars 2 and 3 of interest if( (sv == 2).or.(sv == 3) ) then ! only state vars 2 and 3 of interest
read (unit,610,END=680) line ! read line with value of state var read (unit,610,END=620) line ! read line with value of state var
pos = IO_stringPos(line,1) pos = IO_stringPos(line,1)
do while (scan(IO_stringValue(line,pos,1),'+-',back=.true.)>1) ! is noEfloat value? do while (scan(IO_stringValue(line,pos,1),'+-',back=.true.)>1) ! is noEfloat value?
val = NINT(IO_fixedNoEFloatValue (line,(/0,20/),1)) ! state var's value val = NINT(IO_fixedNoEFloatValue (line,(/0,20/),1)) ! state var's value
mesh_maxValStateVar(sv-1) = max(val,mesh_maxValStateVar(sv-1)) ! remember max val of material and texture index mesh_maxValStateVar(sv-1) = max(val,mesh_maxValStateVar(sv-1)) ! remember max val of material and texture index
if (initialcondTableStyle == 2) then if (initialcondTableStyle == 2) then
read (unit,610,END=680) line ! read extra line read (unit,610,END=620) line ! read extra line
read (unit,610,END=680) line ! read extra line read (unit,610,END=620) line ! read extra line
end if endif
contInts = IO_continousIntValues(unit,mesh_Nelems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) ! get affected elements contInts = IO_continousIntValues(unit,mesh_Nelems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) ! get affected elements
do i = 1,contInts(1) do i = 1,contInts(1)
e = mesh_FEasCP('elem',contInts(1+i)) e = mesh_FEasCP('elem',contInts(1+i))
mesh_element(1+sv,e) = val mesh_element(1+sv,e) = val
end do enddo
if (initialcondTableStyle == 0) read (unit,610,END=680) line ! ignore IP range for old table style if (initialcondTableStyle == 0) read (unit,610,END=620) line ! ignore IP range for old table style
read (unit,610,END=680) line read (unit,610,END=620) line
pos = IO_stringPos(line,1) pos = IO_stringPos(line,1)
end do enddo
end if endif
else else
read (unit,610,END=680) line read (unit,610,END=620) line
end if endif
end do enddo
680 return 620 return
END SUBROUTINE END SUBROUTINE
@ -2673,8 +2661,9 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
use IO use IO
implicit none implicit none
integer(pint) unit,i,j,n,e,t,NnodesStillToBeRead integer(pInt), parameter :: maxNchunks = 66
integer(pInt), dimension (133) :: pos integer(pint) unit,i,j,n,e
integer(pInt), dimension (1+2*maxNchunks) :: pos
integer(pInt), dimension (:), allocatable :: node_seen integer(pInt), dimension (:), allocatable :: node_seen
character*300 line character*300 line
@ -2686,47 +2675,32 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
rewind(unit) rewind(unit)
do do
read (unit,610,END=690) line read (unit,610,END=620) line
pos = IO_stringPos(line,1) pos = IO_stringPos(line,1)
if( IO_lc(IO_stringValue(line,pos,1)) == 'connectivity' ) then if( IO_lc(IO_stringValue(line,pos,1)) == 'connectivity' ) then
read (unit,610,END=690) line ! Garbage line read (unit,610,END=620) line ! Garbage line
readThisElement: do i=1,mesh_Nelems ! read all elements do i=1,mesh_Nelems
NnodesStillToBeRead = -1_pInt ! number of nodes per element still to be read read (unit,610,END=620) line
do while (NnodesStillToBeRead /= 0) pos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max (plus ID, type)
read (unit,610,END=690) line e = mesh_FEasCP('elem',IO_IntValue(line,pos,1))
pos = IO_stringPos(line,2+FE_maxNnodes) ! limit to maximal number of nodes (plus ID, type) if (e /= 0) then ! disregard non CP elems
if (NnodesStillToBeRead < 0) then ! first line of element with FE ID and element type node_seen = 0_pInt
e = mesh_FEasCP('elem',IO_intValue(line,pos,1)) do j = 1,FE_Nnodes(FE_mapElemtype(IO_StringValue(line,pos,2)))
if (e < 1) exit readThisElement ! disregard non-CP elems n = mesh_FEasCP('node',IO_IntValue (line,pos,j+2))
t = FE_mapElemtype(IO_StringValue(line,pos,2)) ! element type if (all(node_seen /= n)) then
NnodesStillToBeRead = FE_Nnodes(t) mesh_sharedElem(1,n) = mesh_sharedElem(1,n) + 1
node_seen = 0_pInt mesh_sharedElem(1+mesh_sharedElem(1,n),n) = e
do j = 1,(pos(1)-2) end if
n = mesh_FEasCP('node',IO_IntValue (line,pos,j+2)) node_seen(j) = n
if (all(node_seen /= n)) then enddo
mesh_sharedElem(1,n) = mesh_sharedElem(1,n) + 1 call IO_skipChunks(unit,FE_NoriginalNodes(mesh_element(2,e))-(pos(1)-2)) ! read on if FE_Nnodes exceeds node count present on current line
mesh_sharedElem(1+mesh_sharedElem(1,n),n) = e end if
end if end do
node_seen(j) = n
NnodesStillToBeRead = NnodesStillToBeRead - 1
end do
else
do j = 1,pos(1)
n = mesh_FEasCP('node',IO_IntValue (line,pos,j))
if (all(node_seen /= n)) then
mesh_sharedElem(1,n) = mesh_sharedElem(1,n) + 1
mesh_sharedElem(1+mesh_sharedElem(1,n),n) = e
end if
node_seen(FE_Nnodes(t)-NnodesStillToBeRead+1) = n
NnodesStillToBeRead = NnodesStillToBeRead - 1
end do
end if
end do
end do readThisElement
exit exit
end if end if
end do end do
690 return
620 return
END SUBROUTINE END SUBROUTINE
@ -2827,7 +2801,7 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
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 true subnodes do n = 1,FE_NsubNodes(t) ! now for the true subnodes
do p = 1,FE_Nips(t) ! loop through parents do p = 1,FE_Nips(t) ! loop through possible parent nodes
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) + &
@ -2854,7 +2828,7 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
use math, only: math_volTetrahedron use math, only: math_volTetrahedron
implicit none implicit none
integer(pInt) e,f,t,i,j,k,n,test 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), 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
@ -2863,7 +2837,7 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
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
write(6,'(a10,x,a20,x,a20,x,a25,3(x,a6))') 'FE_Nips','FE_NipNeighbors','FE_NipFaceNodes','FE_subNodeOnIPFace','x','y','z' !write(6,'(a10,x,a20,x,a20,x,a25,3(x,a6))') 'FE_Nips','FE_NipNeighbors','FE_NipFaceNodes','FE_subNodeOnIPFace','x','y','z'
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
@ -2874,16 +2848,16 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
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)
write(6,'(i10,x,i20,x,i20,x,i25,3(x,f6.3))') i,f,n,FE_subNodeOnIPFace(n,f,i,t),& write(6,'(i10,x,i20,x,i20,x,i25,3(x,f6.3))') i,f,n,FE_subNodeOnIPFace(n,f,i,t),&
mesh_subNodeCoord(1,FE_subNodeOnIPFace(n,f,i,t),e),& mesh_subNodeCoord(1,FE_subNodeOnIPFace(n,f,i,t),e),&
mesh_subNodeCoord(2,FE_subNodeOnIPFace(n,f,i,t),e),& mesh_subNodeCoord(2,FE_subNodeOnIPFace(n,f,i,t),e),&
mesh_subNodeCoord(3,FE_subNodeOnIPFace(n,f,i,t),e) mesh_subNodeCoord(3,FE_subNodeOnIPFace(n,f,i,t),e)
end do end do
end do end do
do j = 1,mesh_maxNnodes+mesh_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+1,mesh_maxNnodes+mesh_maxNsubNodes ! walk through remainder of list do k = j+1,mesh_maxNnodes+mesh_maxNsubNodes ! walk through remainder of list
if (all((gravityNodePos(:,j) - gravityNodePos(:,k)) == 0.0_pReal)) then ! found match if (all(abs(gravityNodePos(:,j) - gravityNodePos(:,k)) < 1.0e-100_pReal)) then ! found duplicate
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
@ -3003,7 +2977,7 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
do i = 1,FE_Nips(mesh_element(2,e)) do i = 1,FE_Nips(mesh_element(2,e))
write (6,"(i5,x,i5,x,e15.8)") e,i,mesh_IPvolume(i,e) write (6,"(i5,x,i5,x,e15.8)") e,i,mesh_IPvolume(i,e)
do f = 1,FE_NipNeighbors(mesh_element(2,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) 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 enddo
enddo enddo