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 createtiedsurfs(nodface,ipoface,set,istartset,
20     &     iendset,ialset,tieset,inomat,ne,ipkon,lakon,kon,ntie,
21     &     tietol,nalset,nk,nset,iactive)
22!
23      implicit none
24!
25!     creates ties for the surfaces in between domains of an
26!     electromagnetic calculation
27!
28      character*8 lakon(*)
29      character*81 set(*),tieset(3,*),surfset
30!
31      integer ipkon(*),kon(*),ne,nodface(5,*),ipoface(*),istartset(*),
32     &     iendset(*),ialset(*),inomat(*),ithree,ifour,ifaceq(8,6),iset,
33     &     ifacet(6,4),ifacew(8,5),ifree,ifreenew,index,indexold,id,
34     &     i,j,k,iactive(3),ntie,nodes(4),nalset,nk,nset,indexe
35!
36      real*8 tietol(3,*)
37!
38!     nodes belonging to the element faces
39!
40      data ifaceq /4,3,2,1,11,10,9,12,
41     &     5,6,7,8,13,14,15,16,
42     &     1,2,6,5,9,18,13,17,
43     &     2,3,7,6,10,19,14,18,
44     &     3,4,8,7,11,20,15,19,
45     &     4,1,5,8,12,17,16,20/
46      data ifacet /1,3,2,7,6,5,
47     &     1,2,4,5,9,8,
48     &     2,3,4,6,10,9,
49     &     1,4,3,8,10,7/
50      data ifacew /1,3,2,9,8,7,0,0,
51     &     4,5,6,10,11,12,0,0,
52     &     1,2,5,4,7,14,10,13,
53     &     2,3,6,5,8,15,11,14,
54     &     4,6,3,1,12,15,9,13/
55!
56      ithree=3
57      ifour=4
58!
59!     determining the external element faces of the electromagnetic mesh
60!     the faces are catalogued by the three lowes nodes numbers
61!     in ascending order. ipoface(i) points to a face for which
62!     node i is the lowest node and nodface(1,ipoface(i)) and
63!     nodface(2,ipoface(i)) are the next lower ones.
64!     nodface(3,ipoface(i)) contains the element number,
65!     nodface(4,ipoface(i)) the face number and nodface(5,ipoface(i))
66!     is a pointer to the next surface for which node i is the
67!     lowest node; if there are no more such surfaces the pointer
68!     has the value zero
69!     An external element face is one which belongs to one element
70!     only
71!
72      ifree=1
73      do i=1,6*ne-1
74        nodface(5,i)=i+1
75      enddo
76      do i=1,ne
77        if(ipkon(i).lt.0) cycle
78        if(lakon(i)(1:3).ne.'C3D') cycle
79        indexe=ipkon(i)
80        if((lakon(i)(4:4).eq.'2').or.(lakon(i)(4:4).eq.'8')) then
81          do j=1,6
82            do k=1,4
83              nodes(k)=kon(indexe+ifaceq(k,j))
84            enddo
85            call insertsorti(nodes,ifour)
86c     call isortii(nodes,iaux,ifour,kflag)
87            indexold=0
88            index=ipoface(nodes(1))
89            do
90!
91!     adding a surface which has not been
92!     catalogued so far
93!
94              if(index.eq.0) then
95                ifreenew=nodface(5,ifree)
96                nodface(1,ifree)=nodes(2)
97                nodface(2,ifree)=nodes(3)
98                nodface(3,ifree)=i
99                nodface(4,ifree)=j
100                nodface(5,ifree)=ipoface(nodes(1))
101                ipoface(nodes(1))=ifree
102                ifree=ifreenew
103                exit
104              endif
105!
106!     removing a surface which has already
107!     been catalogued
108!
109              if((nodface(1,index).eq.nodes(2)).and.
110     &             (nodface(2,index).eq.nodes(3))) then
111                if(indexold.eq.0) then
112                  ipoface(nodes(1))=nodface(5,index)
113                else
114                  nodface(5,indexold)=nodface(5,index)
115                endif
116                nodface(5,index)=ifree
117                ifree=index
118                exit
119              endif
120              indexold=index
121              index=nodface(5,index)
122            enddo
123          enddo
124        elseif((lakon(i)(4:4).eq.'4').or.(lakon(i)(4:5).eq.'10')) then
125          do j=1,4
126            do k=1,3
127              nodes(k)=kon(indexe+ifacet(k,j))
128            enddo
129            call insertsorti(nodes,ithree)
130c     call isortii(nodes,iaux,ithree,kflag)
131            indexold=0
132            index=ipoface(nodes(1))
133            do
134!
135!     adding a surface which has not been
136!     catalogued so far
137!
138              if(index.eq.0) then
139                ifreenew=nodface(5,ifree)
140                nodface(1,ifree)=nodes(2)
141                nodface(2,ifree)=nodes(3)
142                nodface(3,ifree)=i
143                nodface(4,ifree)=j
144                nodface(5,ifree)=ipoface(nodes(1))
145                ipoface(nodes(1))=ifree
146                ifree=ifreenew
147                exit
148              endif
149!
150!     removing a surface which has already
151!     been catalogued
152!
153              if((nodface(1,index).eq.nodes(2)).and.
154     &             (nodface(2,index).eq.nodes(3))) then
155                if(indexold.eq.0) then
156                  ipoface(nodes(1))=nodface(5,index)
157                else
158                  nodface(5,indexold)=nodface(5,index)
159                endif
160                nodface(5,index)=ifree
161                ifree=index
162                exit
163              endif
164              indexold=index
165              index=nodface(5,index)
166            enddo
167          enddo
168        else
169          do j=1,5
170            if(j.le.2) then
171              do k=1,3
172                nodes(k)=kon(indexe+ifacew(k,j))
173              enddo
174              call insertsorti(nodes,ithree)
175c     call isortii(nodes,iaux,ithree,kflag)
176            else
177              do k=1,4
178                nodes(k)=kon(indexe+ifacew(k,j))
179              enddo
180              call insertsorti(nodes,ifour)
181c     call isortii(nodes,iaux,ifour,kflag)
182            endif
183            indexold=0
184            index=ipoface(nodes(1))
185            do
186!
187!     adding a surface which has not been
188!     catalogued so far
189!
190              if(index.eq.0) then
191                ifreenew=nodface(5,ifree)
192                nodface(1,ifree)=nodes(2)
193                nodface(2,ifree)=nodes(3)
194                nodface(3,ifree)=i
195                nodface(4,ifree)=j
196                nodface(5,ifree)=ipoface(nodes(1))
197                ipoface(nodes(1))=ifree
198                ifree=ifreenew
199                exit
200              endif
201!
202!     removing a surface which has already
203!     been catalogued
204!
205              if((nodface(1,index).eq.nodes(2)).and.
206     &             (nodface(2,index).eq.nodes(3))) then
207                if(indexold.eq.0) then
208                  ipoface(nodes(1))=nodface(5,index)
209                else
210                  nodface(5,indexold)=nodface(5,index)
211                endif
212                nodface(5,index)=ifree
213                ifree=index
214                exit
215              endif
216              indexold=index
217              index=nodface(5,index)
218            enddo
219          enddo
220        endif
221      enddo
222!
223!     boundary of domain 1 (phi-domain)
224!
225      surfset(1:22)='ELECTROMAGNETICSZONE1T'
226      do i=23,81
227        surfset(i:i)=' '
228      enddo
229      call cident81(set,surfset,nset,id)
230      nset=nset+1
231      do j=nset,id+2,-1
232        istartset(j)=istartset(j-1)
233        iendset(j)=iendset(j-1)
234        set(j)=set(j-1)
235      enddo
236      iset=id+1
237ccc   to remove start
238c     iset=nset
239ccc   to remove end
240      istartset(iset)=nalset+1
241      do i=1,nk
242        if(ipoface(i).eq.0) cycle
243        if(inomat(i).ne.1) cycle
244        index=ipoface(i)
245        do
246          if(index.eq.0) exit
247          nalset=nalset+1
248          ialset(nalset)=10*nodface(3,index)+nodface(4,index)
249          index=nodface(5,index)
250        enddo
251      enddo
252      if(istartset(iset).gt.nalset) then
253        iactive(1)=0
254        do j=iset+1,nset
255          istartset(j-1)=istartset(j)
256          iendset(j-1)=iendset(j)
257          set(j-1)=set(j)
258        enddo
259        nset=nset-1
260      else
261        iactive(1)=iset
262ccc   to remove start
263c     set(iset)(1:22)='ELECTROMAGNETICSZONE1T'
264c     do i=23,81
265c     set(iset)(i:i)=' '
266c     enddo
267ccc   to remove end
268        set(iset)=surfset
269        iendset(iset)=nalset
270      endif
271!
272!     boundary of domain 2 (A-V-domain)
273!
274      surfset(1:22)='ELECTROMAGNETICSZONE2T'
275      do i=23,81
276        surfset(i:i)=' '
277      enddo
278      call cident81(set,surfset,nset,id)
279      nset=nset+1
280      do j=nset,id+2,-1
281        istartset(j)=istartset(j-1)
282        iendset(j)=iendset(j-1)
283        set(j)=set(j-1)
284      enddo
285      iset=id+1
286ccc   to remove start
287c     iset=nset
288ccc   to remove end
289      istartset(iset)=nalset+1
290      do i=1,nk
291        if(ipoface(i).eq.0) cycle
292        if(inomat(i).ne.2) cycle
293        index=ipoface(i)
294        do
295          if(index.eq.0) exit
296          nalset=nalset+1
297          ialset(nalset)=10*nodface(3,index)+nodface(4,index)
298          index=nodface(5,index)
299        enddo
300      enddo
301      if(istartset(iset).gt.nalset) then
302        iactive(2)=0
303        do j=iset+1,nset
304          istartset(j-1)=istartset(j)
305          iendset(j-1)=iendset(j)
306          set(j-1)=set(j)
307        enddo
308        nset=nset-1
309      else
310        iactive(2)=iset
311ccc   to remove start
312c     set(iset)(1:22)='ELECTROMAGNETICSZONE2T'
313c     do i=23,81
314c     set(iset)(i:i)=' '
315c     enddo
316ccc   to remove end
317        set(iset)=surfset
318        iendset(iset)=nalset
319      endif
320!
321!     boundary of domain 3 (A-domain)
322!
323      surfset(1:22)='ELECTROMAGNETICSZONE3T'
324      do i=23,81
325        surfset(i:i)=' '
326      enddo
327      call cident81(set,surfset,nset,id)
328      nset=nset+1
329      do j=nset,id+2,-1
330        istartset(j)=istartset(j-1)
331        iendset(j)=iendset(j-1)
332        set(j)=set(j-1)
333      enddo
334      iset=id+1
335ccc   to remove start
336c     iset=nset
337ccc   to remove end
338      istartset(iset)=nalset+1
339      do i=1,nk
340        if(ipoface(i).eq.0) cycle
341        if(inomat(i).ne.3) cycle
342        index=ipoface(i)
343        do
344          if(index.eq.0) exit
345          nalset=nalset+1
346          ialset(nalset)=10*nodface(3,index)+nodface(4,index)
347          index=nodface(5,index)
348        enddo
349      enddo
350      if(istartset(iset).gt.nalset) then
351        iactive(3)=0
352        do j=iset+1,nset
353          istartset(j-1)=istartset(j)
354          iendset(j-1)=iendset(j)
355          set(j-1)=set(j)
356        enddo
357        nset=nset-1
358      else
359        iactive(3)=iset
360ccc   to remove start
361c     set(iset)(1:22)='ELECTROMAGNETICSZONE3T'
362c     do i=23,81
363c     set(iset)(i:i)=' '
364c     enddo
365ccc   to remove end
366        set(iset)=surfset
367        iendset(iset)=nalset
368      endif
369!
370!     create contact ties between the domains
371!
372      do i=1,3
373        do j=1,3
374          if((i.eq.j).or.((i.eq.3).and.(j.eq.2))) cycle
375          if(iactive(i)*iactive(j).gt.0) then
376            ntie=ntie+1
377            tietol(1,ntie)=0.d0
378            tietol(2,ntie)=0.d0
379            write(tieset(1,ntie)(1:1),'(i1)')i
380            write(tieset(1,ntie)(2:2),'(i1)')j
381            do k=3,80
382              tieset(1,ntie)(k:k)=' '
383            enddo
384            tieset(1,ntie)(81:81)='E'
385!
386!     slave set
387!     this set does not have to be identified as nodal
388!     or facial at this stage
389!
390            tieset(2,ntie)(1:20)='ELECTROMAGNETICSZONE'
391            write(tieset(2,ntie)(21:21),'(i1)')i
392            do k=22,81
393              tieset(2,ntie)(k:k)=' '
394            enddo
395!
396!     master set
397!
398            tieset(3,ntie)(1:22)='ELECTROMAGNETICSZONE T'
399            write(tieset(3,ntie)(21:21),'(i1)')j
400            do k=23,81
401              tieset(3,ntie)(k:k)=' '
402            enddo
403          endif
404        enddo
405      enddo
406!
407      return
408      end
409