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