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