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