1!
2!     CalculiX - A 3-dimensional finite element program
3!              Copyright (C) 1998-2021 Guido Dhondt
4!
5!     This program is free software; you can redistribute it and/or
6!     modify it under the terms of the GNU General Public License as
7!     published by the Free Software Foundation(version 2);
8!
9!
10!     This program is distributed in the hope that it will be useful,
11!     but WITHOUT ANY WARRANTY; without even the implied warranty of
12!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13!     GNU General Public License for more details.
14!
15!     You should have received a copy of the GNU General Public License
16!     along with this program; if not, write to the Free Software
17!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
18!
19      subroutine findsurface(ipoface,nodface,ne,ipkon,kon,lakon,ntie,
20     &    tieset)
21!
22!     determining the external faces of the mesh and storing
23!     them in fields ipoface and nodface
24!
25      implicit none
26!
27      character*8 lakon(*)
28      character*81 tieset(3,*),slavset
29!
30      integer ipoface(*),nodface(5,*),nodes(4),
31     &  ne,ipkon(*),kon(*),indexe,ifaceq(8,6),ifacet(6,4),index1,
32     &  ifacew(8,5),ithree,ifour,i,j,k,m,
33     &  ifree,index1old,ifreenew,ntie,ipos
34!
35!     nodes belonging to the element faces
36!
37      data ifaceq /4,3,2,1,11,10,9,12,
38     &            5,6,7,8,13,14,15,16,
39     &            1,2,6,5,9,18,13,17,
40     &            2,3,7,6,10,19,14,18,
41     &            3,4,8,7,11,20,15,19,
42     &            4,1,5,8,12,17,16,20/
43      data ifacet /1,3,2,7,6,5,
44     &             1,2,4,5,9,8,
45     &             2,3,4,6,10,9,
46     &             1,4,3,8,10,7/
47      data ifacew /1,3,2,9,8,7,0,0,
48     &             4,5,6,10,11,12,0,0,
49     &             1,2,5,4,7,14,10,13,
50     &             2,3,6,5,8,15,11,14,
51     &             4,6,3,1,12,15,9,13/
52!
53      do m=1,ntie
54!
55!        check for contact conditions
56!
57         if((tieset(1,m)(81:81).eq.'C').or.
58     &      (tieset(1,m)(81:81).eq.'-')) then
59            slavset=tieset(2,m)
60!
61!           check whether facial slave surface;
62!
63            ipos=index(slavset,' ')-1
64            if(slavset(ipos:ipos).eq.'S') then
65               ithree=3
66               ifour=4
67!
68!     determining the external element faces of the solid mesh;
69!     the faces are catalogued by the four lowest node numbers
70!     in ascending order.
71!
72!     ipoface(i) points to a face for which
73!     node i is the lowest end node and nodface(1,ipoface(i)),
74!     nodface(2,ipoface(i)) and nodface(3,ipoface(i)) are the next
75!     lower ones. If the face is triangular nodface(3,ipoface(i))
76!     is zero.
77!
78!     nodface(4,ipoface(i)) contains the face number
79!     (10*element number + local face number) and nodface(5,ipoface(i))
80!     is a pointer to the next surface for which node i is the
81!     lowest node; if there are no more such surfaces the pointer
82!     has the value zero
83!
84!     An external element face is one which belongs to one element
85!     only
86!
87               ifree=1
88               do i=1,6*ne-1
89                  nodface(5,i)=i+1
90               enddo
91               do i=1,ne
92                  if(ipkon(i).lt.0) cycle
93                  if(lakon(i)(1:1).ne.'C') cycle
94                  indexe=ipkon(i)
95!
96!                 hexahedral element
97!
98                  if((lakon(i)(4:4).eq.'2').or.
99     &               (lakon(i)(4:4).eq.'8')) then
100                     do j=1,6
101                        do k=1,4
102                           nodes(k)=kon(indexe+ifaceq(k,j))
103                        enddo
104                        call insertsorti(nodes,ifour)
105c                        call isortii(nodes,iaux,ifour,kflag)
106                        index1old=0
107                        index1=ipoface(nodes(1))
108                        do
109!
110!     adding a surface which has not been
111!     catalogued so far
112!
113                           if(index1.eq.0) then
114                              ifreenew=nodface(5,ifree)
115                              nodface(1,ifree)=nodes(2)
116                              nodface(2,ifree)=nodes(3)
117                              nodface(3,ifree)=nodes(4)
118                              nodface(4,ifree)=10*i+j
119                              nodface(5,ifree)=ipoface(nodes(1))
120                              ipoface(nodes(1))=ifree
121                              ifree=ifreenew
122                              exit
123                           endif
124!
125!     removing a surface which has already
126!     been catalogued
127!
128                           if((nodface(1,index1).eq.nodes(2)).and.
129     &                        (nodface(2,index1).eq.nodes(3)).and.
130     &                        (nodface(3,index1).eq.nodes(4))) then
131                              if(index1old.eq.0) then
132                                 ipoface(nodes(1))=nodface(5,index1)
133                              else
134                                 nodface(5,index1old)=nodface(5,index1)
135                              endif
136                              nodface(5,index1)=ifree
137                              ifree=index1
138                              exit
139                           endif
140                           index1old=index1
141                           index1=nodface(5,index1)
142                        enddo
143                     enddo
144!
145!                 tetrahedral element
146!
147                  elseif((lakon(i)(4:4).eq.'4').or.
148     &                   (lakon(i)(4:5).eq.'10')) then
149                     do j=1,4
150                        do k=1,3
151                           nodes(k)=kon(indexe+ifacet(k,j))
152                        enddo
153                        call insertsorti(nodes,ithree)
154c                        call isortii(nodes,iaux,ithree,kflag)
155                        nodes(4)=0
156                        index1old=0
157                        index1=ipoface(nodes(1))
158                        do
159!
160!     adding a surface which has not been
161!     catalogues so far
162!
163                           if(index1.eq.0) then
164                              ifreenew=nodface(5,ifree)
165                              nodface(1,ifree)=nodes(2)
166                              nodface(2,ifree)=nodes(3)
167                              nodface(3,ifree)=nodes(4)
168                              nodface(4,ifree)=10*i+j
169                              nodface(5,ifree)=ipoface(nodes(1))
170                              ipoface(nodes(1))=ifree
171                              ifree=ifreenew
172                              exit
173                           endif
174!
175!     removing a surface which has already
176!     been catalogued
177!
178                           if((nodface(1,index1).eq.nodes(2)).and.
179     &                        (nodface(2,index1).eq.nodes(3)).and.
180     &                        (nodface(3,index1).eq.nodes(4))) then
181                              if(index1old.eq.0) then
182                                 ipoface(nodes(1))=nodface(5,index1)
183                              else
184                                 nodface(5,index1old)=nodface(5,index1)
185                              endif
186                              nodface(5,index1)=ifree
187                              ifree=index1
188                              exit
189                           endif
190                           index1old=index1
191                           index1=nodface(5,index1)
192                        enddo
193                     enddo
194                  else
195!
196!                 wedge element
197!
198                     do j=1,5
199                        if(j.le.2) then
200                           do k=1,3
201                              nodes(k)=kon(indexe+ifacew(k,j))
202                           enddo
203                           call insertsorti(nodes,ithree)
204c                           call isortii(nodes,iaux,ithree,kflag)
205                           nodes(4)=0
206                        else
207                           do k=1,4
208                              nodes(k)=kon(indexe+ifacew(k,j))
209                           enddo
210                           call insertsorti(nodes,ifour)
211c                           call isortii(nodes,iaux,ifour,kflag)
212                        endif
213                        index1old=0
214                        index1=ipoface(nodes(1))
215                        do
216!
217!     adding a surface which has not been
218!     catalogues so far
219!
220                           if(index1.eq.0) then
221                              ifreenew=nodface(5,ifree)
222                              nodface(1,ifree)=nodes(2)
223                              nodface(2,ifree)=nodes(3)
224                              nodface(3,ifree)=nodes(4)
225                              nodface(4,ifree)=10*i+j
226                              nodface(5,ifree)=ipoface(nodes(1))
227                              ipoface(nodes(1))=ifree
228                              ifree=ifreenew
229                              exit
230                           endif
231!
232!     removing a surface which has already
233!     been catalogued
234!
235                           if((nodface(1,index1).eq.nodes(2)).and.
236     &                        (nodface(2,index1).eq.nodes(3)).and.
237     &                        (nodface(3,index1).eq.nodes(4))) then
238                              if(index1old.eq.0) then
239                                 ipoface(nodes(1))=nodface(5,index1)
240                              else
241                                 nodface(5,index1old)=nodface(5,index1)
242                              endif
243                              nodface(5,index1)=ifree
244                              ifree=index1
245                              exit
246                           endif
247                           index1old=index1
248                           index1=nodface(5,index1)
249                        enddo
250                     enddo
251                  endif
252               enddo
253               exit
254            endif
255         endif
256      enddo
257!
258      return
259      end
260