CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
sutherland_hodgman.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine sutherland_hodgman (nopes, xn, xl2sp, xl2mp, nodem, ipe, ime, iactiveline, nactiveline, ifreeintersec, nelemm, nnodelem, nvertex, pvertex)
 

Function/Subroutine Documentation

◆ sutherland_hodgman()

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