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 real*8 function fuvertex(n,x,cotet,kontet,ipoeln,ieln,node,iedge, 20 & ipoeled,ieled,iedgmid,iedtet) 21! 22! determination of the quality of the ball of linear tetrahedron 23! elements around a node (P.L. George) 24! 25 implicit none 26! 27 integer n,kontet(4,*),ipoeln(*),ieln(2,*),node,ine(2,6),indexe, 28 & ielem,nodes(4),n1,n2,i,j,iedge,ipoeled(*),ieled(2,*), 29 & iedgmid(*),iedtet(6,*) 30! 31 real*8 cotet(3,*),alpha,x(n),volume,surface(4),totsurface, 32 & edgelength(6),radius,hmax,quality 33! 34 data ine /1,2,2,3,1,3,1,4,2,4,3,4/ 35! 36 fuvertex=0.d0 37! 38! alpha is the proporionality factor 39! 40 alpha=dsqrt(6.d0)/12.d0 41! 42 indexe=ipoeln(node) 43! 44 do 45 if(indexe.eq.0) exit 46! 47! an element belonging to the ball of node 48! 49 ielem=ieln(1,indexe) 50! 51 do i=1,4 52 nodes(i)=kontet(i,ielem) 53 if(nodes(i).eq.node) then 54 do j=1,3 55 cotet(j,node)=x(j) 56 enddo 57 endif 58 enddo 59! 60! calculating the volume of the element 61! 62 call calcvol(nodes(1),nodes(2),nodes(3),nodes(4),cotet, 63 & volume) 64 if(volume.le.0.d0) volume=1.d-30 65! 66! calculating area of each face in the element 67! 68 call calcsurf(nodes(1),nodes(2),nodes(3),cotet, 69 & surface(1)) 70 call calcsurf(nodes(2),nodes(3),nodes(4),cotet, 71 & surface(2)) 72 call calcsurf(nodes(3),nodes(4),nodes(1),cotet, 73 & surface(3)) 74 call calcsurf(nodes(4),nodes(1),nodes(2),cotet, 75 & surface(4)) 76! 77! calculating the total surface 78! 79 totsurface=surface(1)+surface(2)+surface(3)+surface(4) 80! 81! radius of the inscribed sphere 82! 83 radius=3.d0*volume/totsurface 84! 85! length of each edge 86! 87 do j=1,6 88 n1=nodes(ine(1,j)) 89 n2=nodes(ine(2,j)) 90 edgelength(j)=dsqrt((cotet(1,n1)-cotet(1,n2))**2+ 91 & (cotet(2,n1)-cotet(2,n2))**2+ 92 & (cotet(3,n1)-cotet(3,n2))**2) 93 enddo 94! 95! maximum edge length 96! 97 hmax=maxval(edgelength) 98! 99! quality of the element 100! 101 quality=alpha*hmax/radius 102! 103! worst quality of elements treated so far 104! 105 fuvertex=max(fuvertex,quality) 106! 107 indexe=ieln(2,indexe) 108 enddo 109 write(*,100) x(1),x(2),x(3),fuvertex 110 100 format('fuvertex ',4(1x,f15.8)) 111! 112 return 113 end 114