1PROGRAM FindParents
2!----------------------------------------------------------------------------
3! Determines parents for boundary elements in Elmer mesh files.
4!
5! Written by: Mikko Lyly 18 May 2004
6! Modified by: Mark Smith 7 March 2014
7!----------------------------------------------------------------------------
8  IMPLICIT NONE
9
10  TYPE HashEntry_t
11     INTEGER :: Element
12     TYPE(HashEntry_t), POINTER :: Next
13  END TYPE HashEntry_t
14
15  TYPE HashTable_t
16     TYPE(HashEntry_t), POINTER :: Head
17  END TYPE HashTable_t
18
19  TYPE(HashTable_t), ALLOCATABLE :: HashTable(:)
20  TYPE(HashEntry_t), POINTER :: HashPtr, HashPtr1
21
22  LOGICAL :: Found
23  INTEGER :: i, j, a(30), Istat
24  INTEGER :: Candidate, Parent(2)
25  INTEGER :: NumberOfParents, Element, CheckCount, NumberOfNodes, Node
26  INTEGER :: Nnodes, Nbulk, Nboundary, BoundaryCode, ElementDim, BoundaryDim
27  INTEGER, ALLOCATABLE :: ElementCode(:)
28
29!---------------------------------------------------------------------------
30
31! Read header:
32! ------------
33  OPEN(10, FILE='mesh.header')
34  READ(10,*) Nnodes, Nbulk, Nboundary
35  CLOSE(10)
36
37! Prepare the element hash table for nodes (inverse connectivity):
38! ----------------------------------------------------------------
39  OPEN(10, FILE='mesh.elements')
40
41  ISTAT = 0
42  ALLOCATE( HashTable( Nnodes ), ElementCode( Nbulk ), STAT = Istat )
43  IF( Istat /= 0 ) THEN
44     PRINT *,'Memory allocation error. Aborting.'
45     STOP
46  END IF
47
48  DO i = 1,Nnodes
49	NULLIFY( HashTable( i ) % Head )
50  ENDDO
51
52  DO i = 1, Nbulk
53     READ(10,*) A(1:3), A(4:3+MOD(A(3),100))
54
55!    A(1) = element number
56!    A(2) = tag
57!    A(3) = code
58!    A(4:) = node numbers
59
60     ElementCode(i) = A(3)
61     NumberOfNodes = MOD(A(3),100)
62
63     DO j = 1,NumberOfNodes
64        Node = A(3+j)
65
66        HashPtr => HashTable( Node ) % Head
67        Found = .FALSE.
68
69        DO WHILE( ASSOCIATED( HashPtr ) )
70           IF( HashPtr % Element == A(1) ) THEN
71              Found = .TRUE.
72              EXIT
73           END IF
74           HashPtr => HashPtr % Next
75        END DO
76
77        IF( .NOT.Found ) THEN
78           ALLOCATE( HashPtr )
79           HashPtr % Element = i
80           HashPtr % Next => HashTable( Node ) % Head
81           HashTable( Node ) % Head => HashPtr
82        END IF
83
84     END DO
85  END DO
86  CLOSE(10)
87
88
89! Make the mesh.boundary -file with parents:
90! ------------------------------------------
91  OPEN(10, FILE='mesh.boundary')
92  OPEN(11, FILE='mesh.boundary.corrected' )
93
94  DO i = 1, Nboundary
95     READ(10,*) A(1:5), A(6:5+MOD(A(5),100))
96
97!    A(1) = boundaryelement number
98!    A(2) = tag
99!    A(3) = left parent (assumed unknown)
100!    A(4) = right parent (assumed unknown)
101!    A(5) = code
102!    A(6:) = node numbers
103
104     BoundaryCode = A(5)
105     NumberOfNodes = MOD(A(5),100)
106
107     SELECT CASE( INT(BoundaryCode/100) )
108     CASE( 1 )
109        BoundaryDim = 0
110     CASE( 2 )
111        BoundaryDim = 1
112     CASE( 3, 4 )
113        BoundaryDim = 2
114     CASE( 5, 6, 7, 8 )
115        BoundaryDim = 3
116     CASE DEFAULT
117        PRINT *,'Cannot detect dimension for boundary element',A(1)
118        BoundaryDim = 0
119     END SELECT
120
121     NumberOfParents = 0
122     HashPtr1 => HashTable( A(6) ) % Head
123     Parent = 0
124
125     DO WHILE( ASSOCIATED( HashPtr1 ) )
126
127        Candidate = HashPtr1 % Element
128        CheckCount = 0
129
130        Do j = 1,NumberOfNodes
131           Node = A(5+j)
132           HashPtr => HashTable( Node ) % Head
133           DO WHILE( ASSOCIATED( HashPtr ) )
134              IF( HashPtr % Element == Candidate ) &
135                   CheckCount = CheckCount+1
136              HashPtr => HashPtr % Next
137           END DO
138        END DO
139
140        IF( CheckCount == NumberOfNodes ) THEN
141           NumberOfParents = NumberOfParents+1
142
143           IF( NumberOfParents > 2 ) THEN
144              PRINT *,'Confused: Found more than 2 parents'
145              PRINT *,'for boundary element =',A(1)
146           END IF
147
148           SELECT CASE( INT(ElementCode(Candidate)/100) )
149           CASE(1)
150              ElementDim = 0
151           CASE(2)
152              ElementDim = 1
153           CASE( 3, 4 )
154              ElementDim = 2
155           CASE( 5, 6, 7, 8 )
156              ElementDim = 3
157           CASE DEFAULT
158              PRINT *,'Cannot detect dimension for element',Candidate
159              ElementDim = 0
160           END SELECT
161
162           IF( ElementDim /= BoundaryDim+1 ) THEN
163              PRINT *,'Confused: Dimension of the boundary element',A(1)
164              PRINT *,'is incompatible with possible parent',Candidate
165           END IF
166
167           IF( NumberOfParents <= 2 ) Parent( NumberOfParents ) = Candidate
168
169        END IF
170
171        HashPtr1 => HashPtr1 % Next
172     END DO
173
174     WRITE(11,'(100I8)') A(1:2), Parent(1:2), A(5), A(6:5+MOD(A(5),100))
175
176  END DO
177
178! Close, deallocate, and destroy hash tables:
179! -------------------------------------------
180  CLOSE(10)
181  CLOSE(11)
182
183  DEALLOCATE( ElementCode )
184
185  DO i = 1,Nnodes
186     HashPtr => HashTable(i) % Head
187     DO WHILE( ASSOCIATED( HashPtr ) )
188        HashPtr1 => HashPtr % Next
189        DEALLOCATE( HashPtr )
190        HashPtr => HashPtr1
191     END DO
192  END DO
193
194! Done:
195! -----
196
197END PROGRAM FindParents
198