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!
20!     Sutherland-Hodgman-Algo for polygon clipping combined with
21!     active line search
22!
23      subroutine sutherland_hodgman(nopes,xn,xl2sp,xl2mp,
24     &  nodem,ipe,ime,iactiveline,nactiveline,ifreeintersec,
25     &  nelemm,nnodelem,nvertex,pvertex)
26!
27      implicit none
28!
29      logical invert,oldactive,altered,border
30!
31      integer nvertex,nopes,ipe(*),ime(4,*),iactiveline(3,*),
32     &  nactiveline,ifreeintersec,itri,nelemm,i,ii,j,k,nnodelem,id,
33     &  nodem(*),ncvertex,node1,node2,modf,node,indexl,ithree,
34     &  insertl(3),ninsertl
35!
36      real*8 pvertex(3,*),xn(3),xl2sp(3,*),pa(3),pb(3),xinters(3),
37     &  xcp(3),diff,dd,xl2mp(3,*),c_pvertex(3,13),t,cedge(3),xtest(3),
38     &     eplane,area,areax,areay,areaz,p1(3),p2(3),areaface
39!
40!
41!
42      data ithree /3/
43!
44      nvertex=0
45      ninsertl=0
46      border=.false.
47!
48!     Initialize Polygon
49!
50      do j=1,nopes
51         nvertex=nvertex+1
52         do k=1,3
53            pvertex(k,nvertex)=xl2sp(k,j)
54         enddo
55      enddo
56        areaface=0.d0
57        do k=1,nvertex-2
58         p1(1)=pvertex(1,k+1)-pvertex(1,1)
59         p1(2)=pvertex(2,k+1)-pvertex(2,1)
60         p1(3)=pvertex(3,k+1)-pvertex(3,1)
61         p2(1)=pvertex(1,k+2)-pvertex(1,1)
62         p2(2)=pvertex(2,k+2)-pvertex(2,1)
63         p2(3)=pvertex(3,k+2)-pvertex(3,1)
64         areax=((p1(2)*p2(3))-(p2(2)*p1(3)))**2
65         areay=(-(p1(1)*p2(3))+(p2(1)*p1(3)))**2
66         areaz=((p1(1)*p2(2))-(p2(1)*p1(2)))**2
67         areaface=areaface+dsqrt(areax+areay+areaz)/2.d0
68       enddo
69!
70!     loop over clipping edges
71!
72      do i=1,(nnodelem)
73         ncvertex=0
74         altered=.false.
75!
76!     generate clipping plane
77!
78         node1=nodem(modf(nnodelem,i))
79         node2=nodem(modf(nnodelem,i+1))
80         invert=.false.
81         if(node2.lt.node1)then
82            node=node1
83            node1=node2
84            node2=node
85            invert=.true.
86         endif
87         indexl=ipe(node1)
88         do
89            if(ime(1,indexl).eq.node2) exit
90            indexl=ime(4,indexl)
91            if(indexl.eq.0) then
92               write(*,*)
93     &        '*ERROR in SH:line was not properly catalogued'
94               write(*,*) 'itri',itri,'node1',node1,'node2',node2
95               call exit(201)
96            endif
97         enddo
98         do k=1,3
99            cedge(k)=xl2mp(k,modf(nnodelem,i+1))
100     &           -xl2mp(k,modf(nnodelem,i))
101         enddo
102         xcp(1)=xn(2)*cedge(3)-xn(3)*cedge(2)
103         xcp(2)=xn(3)*cedge(1)-xn(1)*cedge(3)
104         xcp(3)=xn(1)*cedge(2)-xn(2)*cedge(1)
105         dd=dsqrt(xcp(1)**2+xcp(2)**2+xcp(3)**2)
106         do k=1,3
107            xcp(k)=xcp(k)/dd
108         enddo
109         t=-eplane(xl2mp(1:3,modf(nnodelem,i)),xcp,0.d0)
110!
111!     inside-outside-test
112!
113         do k=1,3
114            xtest(k)=xl2mp(k,modf(nnodelem,i+2))
115         enddo
116         if (eplane(xtest,xcp,t).gt.0)then
117            t=-t
118            do k=1,3
119               xcp(k)=-xcp(k)
120            enddo
121         endif
122         oldactive=.false.
123         call nidentk(iactiveline,indexl, nactiveline,id,ithree)
124         if(id.gt.0.and.iactiveline(1,id).eq.indexl)then
125            oldactive=.true.
126         endif
127         if(oldactive)then
128            nactiveline=nactiveline-1
129            do ii=id,nactiveline
130               do k=1,3
131                  iactiveline(k,ii)=iactiveline(k,ii+1)
132               enddo
133            enddo
134         endif
135!
136         if(nvertex.lt.3) cycle
137!
138!     loop over polygon vertices
139!
140         do j=0, (nvertex-1)
141            do k=1,3
142               pa(k)=pvertex(k,modf(nvertex,j))
143               pb(k)=pvertex(k,modf(nvertex,j+1))
144            enddo
145            if(eplane(pa,xcp,t).le.1.0d-12) then
146               if(eplane(pb,xcp,t).le.1.0d-12)then
147                  ncvertex=ncvertex+1
148                  do k=1,3
149                     c_pvertex(k,ncvertex)=pb(k)
150                  enddo
151               else
152                  if(abs(eplane(pa,xcp,t)).gt.1.0d-10)then
153                     call intersectionpoint(pa,pb,xcp,t,xinters)
154                     diff=(xinters(1)-pa(1))**2+(xinters(2)-pa(2))**2+
155     &                    (xinters(3)-pa(3))**2
156                     diff=dsqrt(diff)
157                     if(diff.gt.1.d-11)then
158                        ncvertex=ncvertex+1
159                        do k=1,3
160                           c_pvertex(k,ncvertex)=xinters(k)
161                        enddo
162                     endif
163                  endif
164!
165                  if((.not.oldactive).and.(.not.altered))then
166                     altered=.true.
167                     nactiveline=nactiveline+1
168                     ifreeintersec=ifreeintersec+1
169                     do ii=nactiveline,id+2,-1
170                        do k=1,3
171                           iactiveline(k,ii)=iactiveline(k,ii-1)
172                        enddo
173                     enddo
174                     iactiveline(1,id+1)=indexl
175                     iactiveline(2,id+1)=nelemm
176                     iactiveline(3,id+1)=ifreeintersec
177                     ninsertl=ninsertl+1
178                     insertl(ninsertl)=indexl
179                  endif
180               endif
181            else
182               if(eplane(pb,xcp,t).le.1.0d-12)then
183                  if(abs(eplane(pb,xcp,t)).lt.1.0d-10)then
184                     do ii=1,3
185                        xinters(ii)=pb(ii)
186                     enddo
187                     ncvertex=ncvertex+1
188                     do k=1,3
189                        c_pvertex(k,ncvertex)=pb(k)
190                     enddo
191                  else
192                     call intersectionpoint(pa,pb,xcp,t,xinters)
193                     ncvertex=ncvertex+2
194                     do k=1,3
195                        c_pvertex(k,ncvertex-1)=xinters(k)
196                        c_pvertex(k,ncvertex)=pb(k)
197                     enddo
198                  endif
199                  if((.not.oldactive).and.(.not.altered))then
200                     if(eplane(pb,xcp,t).lt.0.d0.and.nvertex.gt.2)then
201                        altered=.true.
202                        nactiveline=nactiveline+1
203                        ifreeintersec=ifreeintersec+1
204                        do ii=nactiveline,id+2,-1
205                           do k=1,3
206                              iactiveline(k,ii)=iactiveline(k,ii-1)
207                           enddo
208                        enddo
209                        iactiveline(1,id+1)=indexl
210                        iactiveline(2,id+1)=nelemm
211                        iactiveline(3,id+1)=ifreeintersec
212                        ninsertl=ninsertl+1
213                        insertl(ninsertl)=indexl
214                     endif
215                  endif
216               else
217                  if((.not.oldactive).and.(.not.altered))then
218                     altered=.true.
219                     nactiveline=nactiveline+1
220                     ifreeintersec=ifreeintersec+1
221                     do ii=nactiveline,id+2,-1
222                        do k=1,3
223                           iactiveline(k,ii)=iactiveline(k,ii-1)
224                        enddo
225                     enddo
226                     iactiveline(1,id+1)=indexl
227                     iactiveline(2,id+1)=nelemm
228                     iactiveline(3,id+1)=ifreeintersec
229                     ninsertl=ninsertl+1
230                     insertl(ninsertl)=indexl
231                  endif
232               endif
233            endif
234!
235!     end loop over polygon vertices
236!
237         enddo
238         do j=1,ncvertex
239            do k=1,3
240               pvertex(k,j)=c_pvertex(k,j)
241            enddo
242         enddo
243         nvertex=ncvertex
244!
245!     end loop over clipping edges
246!
247      enddo
248!
249!     remove inserted active lines,if polygon is degenerated
250!
251      if(nvertex.ge.3)then
252         area=0
253        do k=1,nvertex-2
254         p1(1)=pvertex(1,k+1)-pvertex(1,1)
255         p1(2)=pvertex(2,k+1)-pvertex(2,1)
256         p1(3)=pvertex(3,k+1)-pvertex(3,1)
257         p2(1)=pvertex(1,k+2)-pvertex(1,1)
258         p2(2)=pvertex(2,k+2)-pvertex(2,1)
259         p2(3)=pvertex(3,k+2)-pvertex(3,1)
260         areax=((p1(2)*p2(3))-(p2(2)*p1(3)))**2
261         areay=(-(p1(1)*p2(3))+(p2(1)*p1(3)))**2
262         areaz=((p1(1)*p2(2))-(p2(1)*p1(2)))**2
263         area=area+dsqrt(areax+areay+areaz)/2.d0
264       enddo
265         if(border)write(20,*)'border reached'
266      endif
267      if(nvertex.lt.3.or.border)then
268         do i=1,ninsertl
269            oldactive=.false.
270            indexl=insertl(i)
271            call nidentk(iactiveline,indexl, nactiveline,id,ithree)
272            if(id.gt.0)then
273               if(iactiveline(1,id).eq.indexl) oldactive=.true.
274            endif
275
276            if(oldactive)then
277               nactiveline=nactiveline-1
278               do ii=id,nactiveline
279                  do k=1,3
280                     iactiveline(k,ii)=iactiveline(k,ii+1)
281                  enddo
282               enddo
283            endif
284         enddo
285      endif
286!
287      return
288      end
289