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 cattet(kontet,netet_,ifac,ne,ipkon,kon,ifatet,ifreetet,
20     &     bc,itetfa,ifreefa,planfa,ipofa,cotet,cg,ipoeln,ieln,ifreeln,
21     &     lakon,kontetor,iquad,istartset,iendset,ialset,set,nset,filab,
22     &     jfix,iparentel,jobnamec)
23!
24!     catalogueing the tetrahedral elements of the mesh
25!
26      implicit none
27!
28      character*8 lakon(*)
29      character*81 set(*),elset
30      character*87 filab(*)
31      character*132 fnrfn,jobnamec
32!
33      integer kontet(4,*),netet_,i,ifac(4,*),ne,ipkon(*),kon(*),index,j,
34     &     nodes(4),ifatet(4,*),ifreetet,itetfa(2,*),ifreefa,ipofa(*),
35     &     ipoeln(*),ieln(2,*),node,ifreeln,kontetor(6,*),iquad,nset,
36     &     istartset(*),iendset(*),ialset(*),indexe,k,nope,jfix(*),
37     &     iparentel(*),id
38!
39      real*8 bc(4,*),planfa(4,*),cotet(3,*),cg(3,*)
40!
41!     open a file for nodes which are not properly projected on the
42!     free surface
43!
44      open(40,file='WarnNodeNotProjected.nam',status='unknown')
45      write(40,*) '*NSET,NSET=WarnNodeNotProjected'
46!
47!     read the element set name to refine, if any
48!     default: all tetrahedral elements are refined
49!
50      read(filab(48)(27:87),'(a61)') elset(1:61)
51!
52      do i=62,81
53        elset(i:i)=' '
54      enddo
55!
56!     determining the number of the set
57!
58c      do i=1,nset
59c        if(set(i).eq.elset) exit
60c      enddo
61      call cident81(set,elset,nset,id)
62      i=nset+1
63      if(id.gt.0) then
64        if(elset.eq.set(id)) then
65          i=id
66        endif
67      endif
68!
69!     if element set detected:
70!
71      if(i.le.nset) then
72!
73!       replace 'C' by 'A' in all tet elements
74!
75        do j=1,ne
76          if((lakon(j)(1:4).eq.'C3D4').or.
77     &         (lakon(j)(1:5).eq.'C3D10')) then
78            lakon(j)(1:1)='A'
79          endif
80        enddo
81!
82!     replace 'A' by 'C' in all tet elements belonging to the set
83!     to refine
84!
85        do j=istartset(i),iendset(i)
86          if(ialset(j).gt.0) then
87            k=ialset(j)
88            lakon(k)(1:1)='C'
89          else
90            k=ialset(j-2)
91            do
92              k=k-ialset(j)
93              if(k.ge.ialset(j-1)) exit
94              lakon(k)(1:1)='C'
95            enddo
96          endif
97        enddo
98!
99!     setting field jfix in all VERTEX nodes belonging to
100!     elements not to be refined to 1 (including the common
101!     boundaries between the set to be refined and its complement)
102!
103        do i=1,ne
104          if(lakon(i)(1:5).eq.'C3D8I') then
105            nope=8
106          elseif(lakon(i)(4:5).eq.'20') then
107            nope=8
108          elseif(lakon(i)(4:4).eq.'8') then
109            nope=8
110          elseif(lakon(i)(1:5).eq.'C3D10') then
111            cycle
112          elseif(lakon(i)(1:4).eq.'C3D4') then
113            cycle
114          elseif(lakon(i)(4:5).eq.'15') then
115            nope=6
116          elseif(lakon(i)(4:4).eq.'6') then
117            nope=6
118          elseif(lakon(i)(1:2).eq.'ES') then
119            if(lakon(i)(7:7).eq.'C') then
120              cycle
121            else
122              nope=ichar(lakon(i)(8:8))-47
123            endif
124          elseif(lakon(i)(1:4).eq.'MASS') then
125            nope=1
126          elseif(lakon(i)(1:1).eq.'A') then
127            nope=4
128          else
129            cycle
130          endif
131!
132          indexe=ipkon(i)
133          do j=1,nope
134            jfix(kon(indexe+j))=1
135          enddo
136        enddo
137      endif
138!
139!     change on 20200327: original element number should not be
140!     overwritten: needed at reading time for ielmat, ielorien,
141!     ielbody ...
142!
143c!
144c!     determine the first tetrahedral element or first
145c!     unused element
146c!
147c      do i=1,ne
148c        if((lakon(i)(1:4).eq.'C3D4').or.
149c     &       (lakon(i)(1:5).eq.'C3D10').or.
150c     &       (lakon(i)(1:1).eq.char(0))) exit
151c      enddo
152!
153!     determine the first unused element
154!
155      do i=1,ne
156        if(lakon(i)(1:1).eq.char(0)) exit
157      enddo
158c
159      ifreetet=i
160!
161      do
162        do j=i+1,netet_
163!
164!     element number supersedes largest one
165!
166          if(j.gt.ne) then
167            kontet(4,i)=j
168            i=j
169            exit
170          endif
171c!
172c!     tetrahedral element
173c!
174c          if((lakon(j)(1:4).eq.'C3D4').or.
175c     &         (lakon(j)(1:5).eq.'C3D10')) then
176c            kontet(4,i)=j
177c            i=j
178c            exit
179c          endif
180!
181!     unused element
182!
183          if(lakon(j)(1:1).eq.char(0)) then
184            kontet(4,i)=j
185            i=j
186            exit
187          endif
188        enddo
189        if(j.eq.netet_) exit
190      enddo
191      kontet(4,netet_)=0
192!
193!     initialization of ipofa and ifac
194!
195      do i=1,4*netet_
196        ifac(4,i)=i+1
197      enddo
198      ifac(4,4*netet_)=0
199!
200!     initialization of ieln
201!
202      do i=1,4*netet_
203        ieln(2,i)=i+1
204      enddo
205      ieln(2,4*netet_)=0
206!
207!     adding the tetrahedral elements one by one in kontet
208!     tagging the original elements in kon to be removed
209!     by using a model change in input deck jobname.rfn.inp
210!
211!     creating the file name
212!
213      do i=1,132
214        if(ichar(jobnamec(i:i)).eq.0) exit
215      enddo
216      if(i.gt.125) then
217        write(*,*) '*ERROR in cattet'
218        write(*,*) '       jobname has more than 124 characters:'
219        write(*,*) jobnamec(1:132)
220        call exit(201)
221      endif
222      fnrfn(1:i+7)=jobnamec(1:i-1)//'.rfn.inp'
223!
224!     opening the file
225!
226      open(2,file=fnrfn(1:i+7),status='unknown')
227!
228      do i=1,ne
229        if(ipkon(i).lt.0) cycle
230        if((lakon(i)(1:4).ne.'C3D4').and.
231     &       (lakon(i)(1:5).ne.'C3D10')) cycle
232        index=ipkon(i)
233        do j=1,4
234          nodes(j)=kon(index+j)
235        enddo
236        write(2,101)
237        write(2,102) i
238!
239!     if C3D10: store the middle nodes
240!     If there is at least one C3D10 element iquad is set to 1
241!     which means that the refined mesh will be fully quadratic
242!
243        if(lakon(i)(4:4).eq.'1') then
244          iquad=1
245          do j=1,6
246            kontetor(j,ifreetet)=kon(index+4+j)
247          enddo
248        endif
249!
250!       defining the initial parent element (= element itself)
251!
252        iparentel(ifreetet)=i
253        call generatetet_refine(kontet,ifatet,ifreetet,bc,ifac,itetfa,
254     &       ifreefa,planfa,ipofa,nodes,cotet,cg)
255      enddo
256!
257 101  format('*MODEL CHANGE,TYPE=ELEMENT,REMOVE')
258 102  format(i10)
259      close(2)
260!
261!     generating the element per node relationship
262!
263      do j=1,netet_
264        if(kontet(1,j).eq.0) cycle
265        do i=1,4
266          node=kontet(i,j)
267          index=ifreeln
268          ieln(1,index)=j
269          ifreeln=ieln(2,index)
270          if(ifreeln.eq.0) then
271            write(*,*) '*ERROR in cattet: increase the'
272            write(*,*) '       dimension of ieln'
273            call exit(201)
274          endif
275          ieln(2,index)=ipoeln(node)
276          ipoeln(node)=index
277        enddo
278      enddo
279!
280      return
281      end
282
283