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 fumid(n,x,cotet,kontet,ipoeln,ieln,node,iedge, 20 & ipoeled,ieled,iedgmid,iedtet) 21! 22! determination of the quality of the shell of quadratic tetrahedron 23! elements around a node; consists of a combination of: 24! - the quality of the linear tetrahedrons in which it can be 25! subdivided (P.L. George) 26! - a penalty function for negative Jacobians at the integration 27! points 28! 29 implicit none 30! 31 integer n,kontet(4,*),ipoeln(*),ieln(2,*),node,ine(2,6),indexe, 32 & ielem,nodes(10),n1,n2,i,j,k,m,iedge,ipoeled(*),ieled(2,*), 33 & iedgmid(*),iedtet(6,*),i6(4,8),nodesedge(4),iflag 34! 35 real*8 cotet(3,*),alpha,x(n),volume,surface(4),totsurface, 36 & edgelength(6),radius,hmax,quality,xi,et,ze,shp(4,10),xsj, 37 & xl(3,10) 38! 39! nodes belonging to the linear sub-tetrahedra 40! 41 data i6 /8,9,10,4,1,5,7,8,7,6,3,10,9,8,10,7, 42 & 8,9,5,7,9,10,6,7,5,6,7,9,5,2,6,9/ 43! 44! edges of a tetrahedron 45! 46 data ine /1,2,2,3,1,3,1,4,2,4,3,4/ 47! 48 data iflag /2/ 49! 50 include "gauss.f" 51! 52 fumid=0.d0 53! 54! alpha is the proporionality factor 55! 56 alpha=dsqrt(6.d0)/12.d0 57! 58 indexe=ipoeled(iedge) 59! 60 do 61 if(indexe.eq.0) exit 62 ielem=ieled(1,indexe) 63! 64! determine the coordinates of the nodes belonging to the element 65! 66! vertex nodes 67! 68 do i=1,4 69 nodes(i)=kontet(i,ielem) 70 do k=1,3 71 xl(k,i)=cotet(k,nodes(i)) 72 enddo 73 enddo 74! 75! middle nodes 76! 77 do i=1,6 78 nodes(i+4)=iedgmid(iedtet(i,ielem)) 79 do k=1,3 80 xl(k,i+4)=cotet(k,nodes(i+4)) 81 enddo 82! 83! replace the coordinates by the optimization variables 84! 85 if(nodes(i+4).eq.node) then 86 do m=1,3 87 cotet(m,node)=x(m) 88 xl(m,i+4)=x(m) 89 enddo 90 endif 91 enddo 92! 93! loop over the linear sub-tetrahedrons 94! 95 do i=1,8 96! 97 call calcvol(nodes(i6(1,i)),nodes(i6(2,i)),nodes(i6(3,i)), 98 & nodes(i6(4,i)),cotet,volume) 99 if(volume.le.1.d-15) volume=1.d-15 100! 101! calculating the area of each face of the element 102! 103 call calcsurf(nodes(i6(1,i)),nodes(i6(2,i)),nodes(i6(3,i)), 104 & cotet,surface(1)) 105 call calcsurf(nodes(i6(2,i)),nodes(i6(3,i)),nodes(i6(4,i)), 106 & cotet,surface(2)) 107 call calcsurf(nodes(i6(3,i)),nodes(i6(4,i)),nodes(i6(1,i)), 108 & cotet,surface(3)) 109 call calcsurf(nodes(i6(4,i)),nodes(i6(1,i)),nodes(i6(2,i)), 110 & cotet,surface(4)) 111! 112! calculating the total surface 113! 114 totsurface=surface(1)+surface(2)+surface(3)+surface(4) 115! 116! radius of the inscribed sphere 117! 118 radius=3.d0*volume/totsurface 119! 120! length of each edge 121! 122 nodesedge(1)=nodes(i6(1,i)) 123 nodesedge(2)=nodes(i6(2,i)) 124 nodesedge(3)=nodes(i6(3,i)) 125 nodesedge(4)=nodes(i6(4,i)) 126! 127 do j=1,6 128 n1=nodesedge(ine(1,j)) 129 n2=nodesedge(ine(2,j)) 130 edgelength(j)=dsqrt((cotet(1,n1)-cotet(1,n2))**2+ 131 & (cotet(2,n1)-cotet(2,n2))**2+ 132 & (cotet(3,n1)-cotet(3,n2))**2) 133 enddo 134! 135! maximum edge length 136! 137 hmax=maxval(edgelength) 138! 139! quality of the element 140! 141 quality=alpha*hmax/radius 142! 143 fumid=max(fumid,quality) 144 enddo 145! 146! determine the Jacobian in each integration point 147! notice: the penalty due to a negative Jacobian must be 148! bigger than the size of quality (volume >= 1.d-15) 149! 150 do j=1,4 151 xi=gauss3d5(1,j) 152 et=gauss3d5(2,j) 153 ze=gauss3d5(3,j) 154 call shape10tet(xi,et,ze,xl,xsj,shp,iflag) 155 if(xsj.le.0.d0) then 156 fumid=fumid+1.d30*(1.d0-xsj) 157 exit 158 endif 159 enddo 160! 161 indexe=ieled(2,indexe) 162 enddo 163! 164 return 165 end 166