1!
2!     CalculiX - A 3-dimensional finite element program
3!              Copyright (C) 1998 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 createtet(kontet,ifatet,ielement,inodfa,
20     &  ifreefa,planfa,ipofa,nodes,cotet,iparentelement)
21!
22      implicit none
23!
24      integer nodes(4),nodef(3),kontet(4,*),ifatet(4,*),inodfa(4,*),
25     &  ipofa(*),ifreetet,ifreefa,ifree,ig(3,4),j,
26     &  n1,n2,n3,n4,nxa,nxb,nxc,nxd,nya,nyb,nyc,nyd,nza,nzb,nzc,
27     &  nzd,nx1,nx2,ny1,ny2,nz1,nz2,index,j1,j2,j3,i,n,
28     &  node,indexold,ielement,ndx,ndy,ndz,iparentelement
29!
30      real*8 planfa(4,*),cotet(3,*),dd
31!
32      data ig /2,3,4,3,4,1,4,1,2,1,2,3/
33!
34!     updating the node per element relationship
35!
36      do i=1,4
37         kontet(i,ielement)=nodes(i)
38      enddo
39!
40!     creating faces
41!
42      do i=1,4
43         nodef(1)=nodes(ig(1,i))
44         nodef(2)=nodes(ig(2,i))
45         nodef(3)=nodes(ig(3,i))
46!
47         n=3
48         call insertsorti(nodef,n)
49c         call isortii(nodef,idum,n,kflag)
50!
51!        check whether face already exists
52!
53         node=nodef(1)
54         index=ipofa(node)
55!
56         do
57            if(index.eq.0) exit
58            if((inodfa(2,index).eq.nodef(2)).and.
59     &         (inodfa(3,index).eq.nodef(3))) exit
60            indexold=index
61            index=inodfa(4,index)
62         enddo
63!
64         if(index.eq.0) then
65            index=ifreefa
66            ifreefa=inodfa(4,ifreefa)
67            if(ifreefa.eq.0) then
68               write(*,*) '*ERROR in generatet: increase the dimension'
69               write(*,*) '       of inodfa'
70            endif
71            inodfa(1,index)=nodef(1)
72            inodfa(2,index)=nodef(2)
73            inodfa(3,index)=nodef(3)
74            inodfa(4,index)=0
75            if(ipofa(node).eq.0) then
76               ipofa(node)=index
77            else
78               inodfa(4,indexold)=index
79            endif
80!
81            call planeeq(cotet,nodef,planfa(1,index))
82!
83         endif
84!
85!        the face number in ifatet is negative, if the equation
86!        of the face plane is such, that its value in the
87!        remaining node of the tetrahedron is negative
88!
89         dd=planfa(1,index)*cotet(1,nodes(i))+
90     &      planfa(2,index)*cotet(2,nodes(i))+
91     &      planfa(3,index)*cotet(3,nodes(i))+
92     &      planfa(4,index)
93         if(dabs(dd).lt.1.d-10) then
94            write(*,*) '*WARNING in createtet: element ',
95     &                     iparentelement
96            write(*,*) '         is extremely flat'
97            write(*,*) '         the element is deleted'
98            ielement=ielement-1
99            return
100         endif
101         if(dd.ge.0.d0) then
102            ifatet(i,ielement)=index
103         else
104            ifatet(i,ielement)=-index
105         endif
106      enddo
107!
108      return
109      end
110