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 getnodesinitetmesh(ne,lakon,ipkon,kon,istartset,
20     &     iendset,ialset,set,nset,filab,inodestet,nnodestet,
21     &     nodface,ipoface,nk)
22!
23!     reading the initial tet mesh which will be refined
24!
25      implicit none
26!
27      character*8 lakon(*)
28      character*81 set(*),elset
29      character*87 filab(*)
30!
31      integer ipkon(*),kon(*),istartset(*),iendset(*),ialset(*),
32     &     inodestet(*),nnodestet,i,j,k,m,node,ne,nset,indexe,id,
33     &     nodface(5,*),ipoface(*),nodes(3),ifree,ithree,kflag,indexold,
34     &     index,ifreenew,iel,iset,j1,nk,iaux,ifacet(6,4)
35!
36!     nodes belonging to the element faces
37!
38      data ifacet /1,3,2,7,6,5,
39     &     1,2,4,5,9,8,
40     &     2,3,4,6,10,9,
41     &     1,4,3,8,10,7/
42!
43      kflag=1
44      ithree=3
45!
46!     identify the set name/s if they exist
47!
48      read(filab(48)(27:87),'(a61)') elset(1:61)
49!
50      do i=62,81
51        elset(i:i)=' '
52      enddo
53!
54      call cident81(set,elset,nset,id)
55      iset=nset+1
56      if(id.gt.0) then
57        if(elset.eq.set(id)) then
58          iset=id
59        endif
60      endif
61!
62!     determining the external element faces of the unrefined mesh
63!     the faces are catalogued by the three lowes nodes numbers
64!     in ascending order. ipoface(i) points to a face for which
65!     node i is the lowest node and nodface(1,ipoface(i)) and
66!     nodface(2,ipoface(i)) are the next lower ones.
67!     nodface(3,ipoface(i)) contains the element number,
68!     nodface(4,ipoface(i)) the face number and nodface(5,ipoface(i))
69!     is a pointer to the next surface for which node i is the
70!     lowest node; if there are no more such surfaces the pointer
71!     has the value zero
72!     An external element face is one which belongs to one element
73!     only
74!
75      ifree=1
76      do i=1,4*ne-1
77        nodface(5,i)=i+1
78      enddo
79!
80      if(iset.le.nset) then
81        do j1=istartset(iset),iendset(iset)
82          if(ialset(j1).gt.0) then
83            i=ialset(j1)
84!
85!     the elements belonging to the unrefined mesh have been deactivated,
86!     i.e. the transformation ipkon(k)=-2-ipkon(k) was performed in
87!     modelchanges.f
88!
89            indexe=-2-ipkon(i)
90            if(indexe.lt.0) cycle
91            if((lakon(i)(1:4).eq.'C3D4').or.
92     &           (lakon(i)(1:5).eq.'C3D10')) then
93              do j=1,4
94                do k=1,3
95                  nodes(k)=kon(indexe+ifacet(k,j))
96                enddo
97                call isortii(nodes,iaux,ithree,kflag)
98                indexold=0
99                index=ipoface(nodes(1))
100                do
101!
102!     adding a surface which has not been
103!     catalogued so far
104!
105                  if(index.eq.0) then
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            endif
136          endif
137        enddo
138      else
139        do i=1,ne
140!
141!     the elements belonging to the unrefined mesh have been deactivated,
142!     i.e. the transformation ipkon(k)=-2-ipkon(k) was performed in
143!     modelchanges.f
144!
145          indexe=-2-ipkon(i)
146          if(indexe.lt.0) cycle
147          if((lakon(i)(1:4).eq.'C3D4').or.
148     &         (lakon(i)(1:5).eq.'C3D10')) then
149            do j=1,4
150              do k=1,3
151                nodes(k)=kon(indexe+ifacet(k,j))
152              enddo
153              call isortii(nodes,iaux,ithree,kflag)
154              indexold=0
155              index=ipoface(nodes(1))
156              do
157!
158!     adding a surface which has not been
159!     catalogued so far
160!
161                if(index.eq.0) then
162                  ifreenew=nodface(5,ifree)
163                  nodface(1,ifree)=nodes(2)
164                  nodface(2,ifree)=nodes(3)
165                  nodface(3,ifree)=i
166                  nodface(4,ifree)=j
167                  nodface(5,ifree)=ipoface(nodes(1))
168                  ipoface(nodes(1))=ifree
169                  ifree=ifreenew
170                  exit
171                endif
172!
173!     removing a surface which has already
174!     been catalogued
175!
176                if((nodface(1,index).eq.nodes(2)).and.
177     &               (nodface(2,index).eq.nodes(3))) then
178                  if(indexold.eq.0) then
179                    ipoface(nodes(1))=nodface(5,index)
180                  else
181                    nodface(5,indexold)=nodface(5,index)
182                  endif
183                  nodface(5,index)=ifree
184                  ifree=index
185                  exit
186                endif
187                indexold=index
188                index=nodface(5,index)
189              enddo
190            enddo
191          endif
192        enddo
193      endif
194!
195!     storing the nodes belonging to the external faces in inodestet
196!
197      do i=1,nk
198        index=ipoface(i)
199        do
200          if(index.eq.0) exit
201          iel=nodface(3,index)
202          j=nodface(4,index)
203          indexe=-2-ipkon(iel)
204          if(lakon(iel)(4:4).eq.'4') then
205            do k=1,3
206              node=kon(indexe+ifacet(k,j))
207              call nident(inodestet,node,nnodestet,id)
208              if(id.gt.0) then
209                if(inodestet(id).eq.node) cycle
210              endif
211              nnodestet=nnodestet+1
212              do m=nnodestet,id+2,-1
213                inodestet(m)=inodestet(m-1)
214              enddo
215              inodestet(id+1)=node
216            enddo
217          else
218            do k=1,6
219              node=kon(indexe+ifacet(k,j))
220              call nident(inodestet,node,nnodestet,id)
221              if(id.gt.0) then
222                if(inodestet(id).eq.node) cycle
223              endif
224              nnodestet=nnodestet+1
225              do m=nnodestet,id+2,-1
226                inodestet(m)=inodestet(m-1)
227              enddo
228              inodestet(id+1)=node
229            enddo
230          endif
231          index=nodface(5,index)
232        enddo
233      enddo
234!
235      return
236      end
237
238