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