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