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 subroutine checksharp(nexternedg,iedgextfa,cotet,ifacext,isharp) 20! 21! check which edges of the unrefined mesh are sharp 22! 23! the check is done on the angle between the normals on the 24! adjacent faces 25! 26! To this end middle nodes are NOT taken into account, i.e. quadratic 27! faces are reduced to linear faces 28! 29 implicit none 30! 31 integer nexternedg,iedgextfa(2,*),ifacext(6,*),isharp(*),iflag, 32 & imastfa,i,j,k 33! 34 real*8 cotet(3,*),xi,et,xn1(3),xn2(3),dd,xl(3,3),xsj(3),xs(3,7), 35 & shp(7,3) 36! 37! 38! 39 iflag=2 40 xi=1.d0/3.d0 41 et=1.d0/3.d0 42! 43 do i=1,nexternedg 44! 45! first neighboring face 46! 47 imastfa=iedgextfa(1,i) 48 do j=1,3 49 do k=1,3 50 xl(k,j)=cotet(k,ifacext(j,imastfa)) 51 enddo 52 enddo 53 call shape3tri(xi,et,xl,xsj,xs,shp,iflag) 54 dd=dsqrt(xsj(1)*xsj(1)+xsj(2)*xsj(2)+xsj(3)*xsj(3)) 55 do j=1,3 56 xn1(j)=xsj(j)/dd 57 enddo 58! 59! second neighboring face 60! 61 imastfa=iedgextfa(2,i) 62 do j=1,3 63 do k=1,3 64 xl(k,j)=cotet(k,ifacext(j,imastfa)) 65 enddo 66 enddo 67 call shape3tri(xi,et,xl,xsj,xs,shp,iflag) 68 dd=dsqrt(xsj(1)*xsj(1)+xsj(2)*xsj(2)+xsj(3)*xsj(3)) 69 do j=1,3 70 xn2(j)=xsj(j)/dd 71 enddo 72! 73! if the normals are nearly parallel, the edge is no sharp edge 74! "nearly parallel" means that the angle between the vectors 75! is smaller than 0.0464 degrees. 76! 77 if(dabs(xn1(1)*xn2(1)+xn1(2)*xn2(2)+xn1(3)*xn2(3)-1.d0) 78 & .gt.1.d-10) isharp(i)=1 79 enddo 80! 81 return 82 end 83