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 findextsurface(nodface,ipoface,ne,ipkon,lakon,
20     &  kon,konfa,ipkonfa,nk,lakonfa,nsurfs,ifreemax)
21!
22      implicit none
23!
24!     1) catalogues the external surface of the structure: see
25!     comments further down
26!     2) creates topology fields for the external faces:
27!        ipkonfa(1...nsurfs): pointer into konfa
28!        konfa(..): konfa(ipkonfa(i)+1,....ipkonfa(i+1)) contains
29!                   the nodes belonging to face i
30!        lakonfa(1...nsurfs): label (S3, S4, S6 or S8)
31!
32      character*8 lakon(*),lakonfa(*)
33!
34      integer ipkon(*),kon(*),ne,nodface(5,*),ipoface(*),nk,nopem,m,
35     &  ithree,ifour,ifaceq(8,6),konfa(*),ipkonfa(*),nelemm,jfacem,
36     &  ifacet(6,4),ifacew2(8,5),ifree,ifreenew,index,indexold,
37     &  i,j,k,nodes(4),indexe,konl(26),nope,nsurfs,
38     &  ifacew1(4,5),ifreemax
39!
40!
41!
42!     nodes belonging to the element faces
43!
44      data ifaceq /4,3,2,1,11,10,9,12,
45     &            5,6,7,8,13,14,15,16,
46     &            1,2,6,5,9,18,13,17,
47     &            2,3,7,6,10,19,14,18,
48     &            3,4,8,7,11,20,15,19,
49     &            4,1,5,8,12,17,16,20/
50      data ifacet /1,3,2,7,6,5,
51     &             1,2,4,5,9,8,
52     &             2,3,4,6,10,9,
53     &             1,4,3,8,10,7/
54      data ifacew1 /1,3,2,0,
55     &             4,5,6,0,
56     &             1,2,5,4,
57     &             2,3,6,5,
58     &             3,1,4,6/
59      data ifacew2 /1,3,2,9,8,7,0,0,
60     &             4,5,6,10,11,12,0,0,
61     &             1,2,5,4,7,14,10,13,
62     &             2,3,6,5,8,15,11,14,
63     &             3,1,4,6,9,13,12,15/
64!
65      ithree=3
66      ifour=4
67!
68!     determining the external element faces of the mesh;
69!     the faces are catalogued by the three lowest nodes numbers
70!     in ascending order. ipoface(i) points to a face for which
71!     node i is the lowest node and nodface(1,ipoface(i)) and
72!     nodface(2,ipoface(i)) are the next lower ones.
73!     nodface(3,ipoface(i)) contains the element number,
74!     nodface(4,ipoface(i)) the face number and nodface(5,ipoface(i))
75!     is a pointer to the next surface for which node i is the
76!     lowest node; if there are no more such surfaces the pointer
77!     has the value zero
78!     An external element face is one which belongs to one element
79!     only
80!
81      ifree=1
82      ifreemax=1
83      do i=1,6*ne-1
84         nodface(5,i)=i+1
85      enddo
86      do i=1,ne
87         if(ipkon(i).lt.0) cycle
88         if(lakon(i)(1:3).ne.'C3D') cycle
89         indexe=ipkon(i)
90         if((lakon(i)(4:4).eq.'2').or.(lakon(i)(4:4).eq.'8')) then
91            do j=1,6
92               do k=1,4
93                  nodes(k)=kon(indexe+ifaceq(k,j))
94               enddo
95               call insertsorti(nodes,ifour)
96c               call isortii(nodes,iaux,ifour,kflag)
97               indexold=0
98               index=ipoface(nodes(1))
99               do
100!
101!     adding a surface which has not been
102!     catalogued so far
103!
104                  if(index.eq.0) then
105                     ifreemax=max(ifreemax,ifree)
106                     ifreenew=nodface(5,ifree)
107                     nodface(1,ifree)=nodes(2)
108                     nodface(2,ifree)=nodes(3)
109                     nodface(3,ifree)=i
110                     nodface(4,ifree)=j
111                     nodface(5,ifree)=ipoface(nodes(1))
112                     ipoface(nodes(1))=ifree
113                     ifree=ifreenew
114                     exit
115                  endif
116!
117!     removing a surface which has already
118!     been catalogued
119!
120                  if((nodface(1,index).eq.nodes(2)).and.
121     &                 (nodface(2,index).eq.nodes(3))) then
122                     if(indexold.eq.0) then
123                        ipoface(nodes(1))=nodface(5,index)
124                     else
125                        nodface(5,indexold)=nodface(5,index)
126                     endif
127                     nodface(5,index)=ifree
128                     ifree=index
129                     exit
130                  endif
131                  indexold=index
132                  index=nodface(5,index)
133               enddo
134            enddo
135         elseif((lakon(i)(4:4).eq.'4').or.(lakon(i)(4:5).eq.'10')) then
136            do j=1,4
137               do k=1,3
138                  nodes(k)=kon(indexe+ifacet(k,j))
139               enddo
140               call insertsorti(nodes,ithree)
141c               call isortii(nodes,iaux,ithree,kflag)
142               indexold=0
143               index=ipoface(nodes(1))
144               do
145!
146!     adding a surface which has not been
147!     catalogued so far
148!
149                  if(index.eq.0) then
150                     ifreemax=max(ifreemax,ifree)
151                     ifreenew=nodface(5,ifree)
152                     nodface(1,ifree)=nodes(2)
153                     nodface(2,ifree)=nodes(3)
154                     nodface(3,ifree)=i
155                     nodface(4,ifree)=j
156                     nodface(5,ifree)=ipoface(nodes(1))
157                     ipoface(nodes(1))=ifree
158                     ifree=ifreenew
159                     exit
160                  endif
161!
162!     removing a surface which has already
163!     been catalogued
164!
165                  if((nodface(1,index).eq.nodes(2)).and.
166     &                 (nodface(2,index).eq.nodes(3))) then
167                     if(indexold.eq.0) then
168                        ipoface(nodes(1))=nodface(5,index)
169                     else
170                        nodface(5,indexold)=nodface(5,index)
171                     endif
172                     nodface(5,index)=ifree
173                     ifree=index
174                     exit
175                  endif
176                  indexold=index
177                  index=nodface(5,index)
178               enddo
179            enddo
180         else
181            do j=1,5
182               if(j.le.2) then
183                  do k=1,3
184                     nodes(k)=kon(indexe+ifacew2(k,j))
185                  enddo
186                  call insertsorti(nodes,ithree)
187c                  call isortii(nodes,iaux,ithree,kflag)
188               else
189                  do k=1,4
190                     nodes(k)=kon(indexe+ifacew2(k,j))
191                  enddo
192                  call insertsorti(nodes,ifour)
193c                  call isortii(nodes,iaux,ifour,kflag)
194               endif
195               indexold=0
196               index=ipoface(nodes(1))
197               do
198!
199!     adding a surface which has not been
200!     catalogued so far
201!
202                  if(index.eq.0) then
203                     ifreemax=max(ifreemax,ifree)
204                     ifreenew=nodface(5,ifree)
205                     nodface(1,ifree)=nodes(2)
206                     nodface(2,ifree)=nodes(3)
207                     nodface(3,ifree)=i
208                     nodface(4,ifree)=j
209                     nodface(5,ifree)=ipoface(nodes(1))
210                     ipoface(nodes(1))=ifree
211                     ifree=ifreenew
212                     exit
213                  endif
214!
215!     removing a surface which has already
216!     been catalogued
217!
218                  if((nodface(1,index).eq.nodes(2)).and.
219     &                 (nodface(2,index).eq.nodes(3))) then
220                     if(indexold.eq.0) then
221                        ipoface(nodes(1))=nodface(5,index)
222                     else
223                        nodface(5,indexold)=nodface(5,index)
224                     endif
225                     nodface(5,index)=ifree
226                     ifree=index
227                     exit
228                  endif
229                  indexold=index
230                  index=nodface(5,index)
231               enddo
232            enddo
233         endif
234      enddo
235!
236!     Create fields konfa and ipkonfa for calculation of normals of shell
237!     elements
238!
239      nsurfs=0
240      ifree=0
241      do j=1,nk
242!
243         if(ipoface(j).eq.0) cycle
244         indexe=ipoface(j)
245!
246         do
247            nsurfs=nsurfs+1
248            nelemm=nodface(3,indexe)
249            jfacem=nodface(4,indexe)
250!
251!     treatment of hexahedral elements
252!
253            if(lakon(nelemm)(4:4).eq.'8') then
254               nopem=4
255               nope=8
256            elseif(lakon(nelemm)(4:5).eq.'20') then
257               nopem=8
258               nope=20
259!
260!     treatment of tethrahedral elements
261!
262            elseif(lakon(nelemm)(4:5).eq.'10') then
263               nopem=6
264               nope=10
265            elseif(lakon(nelemm)(4:4).eq.'4') then
266               nopem=3
267               nope=4
268!
269!     treatment of wedge faces
270!
271            elseif(lakon(nelemm)(4:4).eq.'6') then
272               nope=6
273               if(jfacem.le.2) then
274                  nopem=3
275               else
276                  nopem=4
277               endif
278            elseif(lakon(nelemm)(4:5).eq.'15') then
279               nope=15
280               if(jfacem.le.2) then
281                  nopem=6
282               else
283                  nopem=8
284               endif
285            endif
286!
287!     actual position of the nodes
288!
289            do k=1,nope
290               konl(k)=kon(ipkon(nelemm)+k)
291            enddo
292!
293!     quadratic quad shell element
294!
295            if((nope.eq.20).or.
296     &           ((nope.eq.15).and.(jfacem.gt.2))) then
297               lakonfa(nsurfs)='S8'
298               ipkonfa(nsurfs)=ifree
299               do m=1,nopem
300                  if(nope.eq.20) then
301                     konfa(ifree+m)=konl(ifaceq(m,jfacem))
302                  else
303                     konfa(ifree+m)=konl(ifacew2(m,jfacem))
304                  endif
305               enddo
306               ifree=ifree+nopem
307!
308!     linear quad shell element
309!
310            elseif((nope.eq.8).or.
311     &              ((nope.eq.6).and.(jfacem.gt.2))) then
312               lakonfa(nsurfs)='S4'
313               ipkonfa(nsurfs)=ifree
314               do m=1,nopem
315                  if(nope.eq.8) then
316                     konfa(ifree+m)=konl(ifaceq(m,jfacem))
317                  else
318                     konfa(ifree+m)=konl(ifacew1(m,jfacem))
319                  endif
320               enddo
321               ifree=ifree+nopem
322!
323!     quadratic tri shell element
324!
325            elseif((nope.eq.10).or.
326     &              ((nope.eq.15).and.(jfacem.le.2))) then
327               lakonfa(nsurfs)='S6'
328               ipkonfa(nsurfs)=ifree
329               do m=1,nopem
330                  if(nope.eq.10) then
331                     konfa(ifree+m)=konl(ifacet(m,jfacem))
332                  else
333                     konfa(ifree+m)=konl(ifacew2(m,jfacem))
334                  endif
335               enddo
336               ifree=ifree+nopem
337!
338!     linear tri shell element
339!
340            elseif((nope.eq.4).or.
341     &              ((nope.eq.6).and.(jfacem.le.2))) then
342               lakonfa(nsurfs)='S3'
343               ipkonfa(nsurfs)=ifree
344               do m=1,nopem
345                  if(nope.eq.4) then
346                     konfa(ifree+m)=konl(ifacet(m,jfacem))
347                  else
348                     konfa(ifree+m)=konl(ifacew1(m,jfacem))
349                  endif
350               enddo
351               ifree=ifree+nopem
352            endif
353            indexe=nodface(5,indexe)
354            if(indexe.eq.0) exit
355         enddo
356      enddo
357!
358      return
359      end
360