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 generatetet_refine(kontet,ifatet,ifreetet,bc,ifac,
20     &  itetfa,ifreefa,planfa,ipofa,nodes,cotet,cg)
21!
22      implicit none
23!
24      integer nodes(4),nodef(3),kontet(4,*),ifatet(4,*),ifac(4,*),
25     &  itetfa(2,*),ipofa(*),ifreetet,ifreefa,ig(3,4),index,i,n,kflag,
26     &  node,ielement
27!
28      real*8 bc(4,*),planfa(4,*),cotet(3,*),cg(3,*),
29     &  p1x,p1y,p1z,p2x,p2y,p2z,p3x,p3y,p3z,p4x,p4y,p4z,a11,
30     &  a12,a13,a21,a22,a23,a31,a32,a33,d1,d2,d3,det,x,y,z,r,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      ielement=ifreetet
37      ifreetet=kontet(4,ielement)
38      if(ifreetet.eq.0) then
39         write(*,*)
40     &        '*ERROR in generatetet_refine increase the dimension'
41         write(*,*) '       of kontet'
42         call exit(201)
43      endif
44!
45      do i=1,4
46         kontet(i,ielement)=nodes(i)
47      enddo
48!
49!     creating faces
50!
51      do i=1,4
52         nodef(1)=nodes(ig(1,i))
53         nodef(2)=nodes(ig(2,i))
54         nodef(3)=nodes(ig(3,i))
55!
56         n=3
57         kflag=1
58         call insertsorti(nodef,n)
59!
60!        check whether face already exists
61!
62         node=nodef(1)
63         index=ipofa(node)
64!
65         do
66            if(index.eq.0) exit
67            if((ifac(2,index).eq.nodef(2)).and.
68     &         (ifac(3,index).eq.nodef(3))) exit
69            index=ifac(4,index)
70         enddo
71!
72         if(index.eq.0) then
73            index=ifreefa
74            ifreefa=ifac(4,ifreefa)
75            if(ifreefa.eq.0) then
76               write(*,*)
77     &           '*ERROR in generatetet_refine increase the dimension'
78               write(*,*) '       of ifac'
79               call exit(201)
80            endif
81            ifac(1,index)=nodef(1)
82            ifac(2,index)=nodef(2)
83            ifac(3,index)=nodef(3)
84            itetfa(1,index)=ielement
85            ifac(4,index)=ipofa(node)
86            ipofa(node)=index
87!
88            call planeeq(cotet,nodef,planfa(1,index))
89!
90         else
91            if(itetfa(1,index).eq.0) then
92               itetfa(1,index)=ielement
93            else
94               itetfa(2,index)=ielement
95            endif
96         endif
97!
98!        the face number in ifatet is negative, if the equation
99!        of the face plane is such, that its value in the
100!        remaining node of the tetrahedron is negative
101!
102         dd=planfa(1,index)*cotet(1,nodes(i))+
103     &      planfa(2,index)*cotet(2,nodes(i))+
104     &      planfa(3,index)*cotet(3,nodes(i))+
105     &      planfa(4,index)
106         if(dd.ge.0.d0) then
107            ifatet(i,ielement)=index
108         else
109            ifatet(i,ielement)=-index
110         endif
111      enddo
112!
113!     finding the center and radius of the circumscribed sphere
114!
115      p1x=cotet(1,nodes(1))
116      p1y=cotet(2,nodes(1))
117      p1z=cotet(3,nodes(1))
118!
119      p2x=cotet(1,nodes(2))
120      p2y=cotet(2,nodes(2))
121      p2z=cotet(3,nodes(2))
122!
123      p3x=cotet(1,nodes(3))
124      p3y=cotet(2,nodes(3))
125      p3z=cotet(3,nodes(3))
126!
127      p4x=cotet(1,nodes(4))
128      p4y=cotet(2,nodes(4))
129      p4z=cotet(3,nodes(4))
130!
131      a11=p1x-p2x
132      a12=p1y-p2y
133      a13=p1z-p2z
134      d1=((p1x+p2x)*a11+(p1y+p2y)*a12+(p1z+p2z)*a13)/2.d0
135!
136      a21=p1x-p3x
137      a22=p1y-p3y
138      a23=p1z-p3z
139      d2=((p1x+p3x)*a21+(p1y+p3y)*a22+(p1z+p3z)*a23)/2.d0
140!
141      a31=p1x-p4x
142      a32=p1y-p4y
143      a33=p1z-p4z
144      d3=((p1x+p4x)*a31+(p1y+p4y)*a32+(p1z+p4z)*a33)/2.d0
145!
146      det=a11*(a22*a33-a32*a23)-a12*(a21*a33-a31*a23)+
147     &     a13*(a21*a32-a31*a22)
148      x=(d1*(a22*a33-a23*a32)-d2*(a12*a33-a32*a13)+
149     &     d3*(a12*a23-a22*a13))/det
150      y=(-d1*(a21*a33-a31*a23)+d2*(a11*a33-a31*a13)-
151     &     d3*(a11*a23-a21*a13))/det
152      z=(d1*(a21*a32-a31*a22)-d2*(a11*a32-a31*a12)+
153     &       d3*(a11*a22-a21*a12))/det
154!
155      r=dsqrt((x-p1x)**2+(y-p1y)**2+(z-p1z)**2)
156!
157      bc(1,ielement)=x
158      bc(2,ielement)=y
159      bc(3,ielement)=z
160      bc(4,ielement)=r
161!
162!     determining the center of gravity of the element
163!
164      cg(1,ielement)=(p1x+p2x+p3x+p4x)/4.d0
165      cg(2,ielement)=(p1y+p2y+p3y+p4y)/4.d0
166      cg(3,ielement)=(p1z+p2z+p3z+p4z)/4.d0
167!
168      return
169      end
170