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!     center of gravity of the projection of the vertices for
21!     visibility purposes
22!     exact integration for one triangle: routine cubtri
23!     if the surfaces are far enough away, one-point integration
24!     is used
25!
26      subroutine calcview(sideload,vold,co,pmid,e1,e2,e3,
27     &     kontri,nloadtr,adview,auview,dist,idist,area,
28     &     ntrit,mi,jqrad,irowrad,nzsrad,sidemean,ntria,ntrib,
29     &     covered,ng)
30!
31      implicit none
32!
33      character*1 covered(ng,ng)
34!
35      character*20 sideload(*)
36!
37      integer ntr,i,j,k,l,mi(*),ntria,ntrib,
38     &     ncovered,kontri(4,*),nloadtr(*),
39     &     i1,j1,istart,iend,jstart,jend,imin,imid,imax,
40     &     k1,kflag,idist(*),ndist,i2,i3,ng,idi,idj,ntri,
41     &     ix,iy,ntrit,limev,ier,nw,idata(1),ncalls,
42     &     irowrad(*),jqrad(*),nzsrad,i0,nzsradv(3)
43!
44      real*8 w(239),vold(0:mi(2),*),co(3,*),xn(3),xxn,xy(ng),
45     &     pmid(3,*),e3(4,*),e1(3,*),e2(3,*),p1(3),p2(3),p3(3),
46     &     x,y,porigin(3),yqmin,yqmax,xqmin,anglemin,distance,
47     &     xxmid,xqmax,dummy,a(3,3),b(3,3),c(3,3),ddd(3),p31(3),
48     &     xq(3),yq(3),ftij,adview(*),auview(*),dint,dir(3),
49     &     dirloc(3),dist(*),area(*),dd,p21(3),sidemean,r(3,3),
50     &     fform,pl(2,3),epsabs,epsrel,abserr,q(3,3),
51     &     rdata(21),factor,argument
52!
53c      real*8 vertex(3,3),vertexl(2,3),unitvec(3,3)
54!
55      external fform
56!
57      nzsradv(3)=nzsrad
58!
59c      ng=160
60      dint=2.d0/ng
61      anglemin=dacos((ng/2.d0-1.d0)/(ng/2.d0))
62!
63!     factor determines when the numerical integration using cubtri
64!     is replaced by a simplified formula using only the center
65!     of gravity of one of the triangles. The integration over the
66!     other triangle is exact (analytical formula, see
67!     "Radiosity: a Programmer's Perspective", by Ian Ashdown, Wiley, 1994)
68!     If the distance between the center of gravity of the triangles
69!     is larger then factor*the projected sqrt(area) of the triangle on the
70!     hemisphere, the simplified formula is taken
71!
72      factor=0.d0
73!
74      do i=ntria,ntrib
75         if(area(i).lt.1.d-20) cycle
76!
77!        pl contains the coordinates of the vertices of triangle i
78!        in the local coordinate system attached to triangle i
79!
80!        This local coordinate system has its origin in the first node
81!        of triangle i (i1 underneath), its local x-axis along the
82!        edge from node 1 to node 2 and its local y-axis such that
83!        e1 x e2 agrees with the usual corkscrew rule
84!
85         i1=kontri(1,i)
86         if(i1.eq.0) cycle
87         i2=kontri(2,i)
88         i3=kontri(3,i)
89         do j=1,3
90            porigin(j)=co(j,i1)+vold(j,i1)
91            p2(j)=co(j,i2)+vold(j,i2)
92            p3(j)=co(j,i3)+vold(j,i3)
93            p21(j)=p2(j)-porigin(j)
94            p31(j)=p3(j)-porigin(j)
95         enddo
96         pl(1,1)=0.d0
97         pl(2,1)=0.d0
98         pl(1,2)=dsqrt(p21(1)**2+p21(2)**2+p21(3)**2)
99         pl(2,2)=0.d0
100         pl(1,3)=p31(1)*e1(1,i)+p31(2)*e1(2,i)+p31(3)*e1(3,i)
101         pl(2,3)=p31(1)*e2(1,i)+p31(2)*e2(2,i)+p31(3)*e2(3,i)
102!
103c         do k=1,3
104c            unitvec(k,1)=e1(k,i)
105c            unitvec(k,2)=e2(k,i)
106c            unitvec(k,3)=e3(k,i)
107c         enddo
108!
109!     checking which triangles face triangle i
110!
111         ndist=0
112         idi=kontri(4,i)
113         do j=1,ntrit
114            if((kontri(1,j).eq.0).or.(area(j).lt.1.d-20).or.
115     &         (i.eq.j)) cycle
116!
117!           if the angle between the connection of the two triangles
118!           and the base plane of the hemisphere is too small (i.e.
119!           triangle j is too close to the horizon) the viewfactor
120!           for this triangle is not taken into account
121!
122            distance=dsqrt((pmid(1,j)-pmid(1,i))**2+
123     &           (pmid(2,j)-pmid(2,i))**2+
124     &           (pmid(3,j)-pmid(3,i))**2)
125            if((pmid(1,j)*e3(1,i)+pmid(2,j)*e3(2,i)+
126     &           pmid(3,j)*e3(3,i)+e3(4,i))/distance.le.anglemin) cycle
127            if((pmid(1,i)*e3(1,j)+pmid(2,i)*e3(2,j)+
128     &           pmid(3,i)*e3(3,j)+e3(4,j))/distance.le.anglemin) cycle
129c            if(pmid(1,j)*e3(1,i)+pmid(2,j)*e3(2,i)+
130c     &           pmid(3,j)*e3(3,i)+e3(4,i).le.sidemean/800.d0) cycle
131c            if(pmid(1,i)*e3(1,j)+pmid(2,i)*e3(2,j)+
132c     &           pmid(3,i)*e3(3,j)+e3(4,j).le.sidemean/800.d0) cycle
133!
134            idj=kontri(4,j)
135            if(sideload(nloadtr(idi))(18:20).ne.
136     &           sideload(nloadtr(idj))(18:20)) cycle
137!
138            ndist=ndist+1
139            dist(ndist)=distance
140c            dist(ndist)=dsqrt((pmid(1,j)-pmid(1,i))**2+
141c     &           (pmid(2,j)-pmid(2,i))**2+
142c     &           (pmid(3,j)-pmid(3,i))**2)
143            idist(ndist)=j
144         enddo
145         if(ndist.eq.0) cycle
146!
147!     ordering the triangles
148!
149         kflag=2
150         call dsort(dist,idist,ndist,kflag)
151!
152!     initializing the coverage matrix
153!
154         do i1=1,ng
155            xy(i1)=((i1-0.5d0)*dint-1.d0)**2
156         enddo
157!
158         ncovered=0
159         do i1=1,ng
160c            x=((i1-0.5d0)*dint-1.d0)**2
161            x=xy(i1)
162            do j1=1,ng
163c               y=((j1-0.5d0)*dint-1.d0)**2
164c               y=xy(j1)
165               if(x+xy(j1).gt.1.d0) then
166                  covered(i1,j1)='T'
167                  ncovered=ncovered+1
168               else
169                  covered(i1,j1)='F'
170               endif
171            enddo
172         enddo
173!
174         do k1=1,ndist
175            j=idist(k1)
176!
177!           determining the 2-D projection of the vertices
178!           of triangle j
179!
180!           r is the connection of cg of triangle i with the
181!           vertices of triangle j
182!
183!           xq and yq are the local coordinates in the local
184!           coordinate system attached of triangle i of the projection
185!           of the vertices of triangle j on the base of the hemisphere
186!           at triangle i
187!
188            do k=1,3
189               r(k,1)=co(k,kontri(1,j))+vold(k,kontri(1,j))-pmid(k,i)
190            enddo
191            ddd(1)=dsqrt(r(1,1)*r(1,1)+r(2,1)*r(2,1)+r(3,1)*r(3,1))
192            do k=1,3
193               r(k,1)=r(k,1)/ddd(1)
194            enddo
195            xq(1)=r(1,1)*e1(1,i)+r(2,1)*e1(2,i)+r(3,1)*e1(3,i)
196            yq(1)=r(1,1)*e2(1,i)+r(2,1)*e2(2,i)+r(3,1)*e2(3,i)
197!
198            do k=1,3
199               r(k,2)=co(k,kontri(2,j))+vold(k,kontri(2,j))-pmid(k,i)
200            enddo
201            ddd(2)=dsqrt(r(1,2)*r(1,2)+r(2,2)*r(2,2)+r(3,2)*r(3,2))
202            do k=1,3
203               r(k,2)=r(k,2)/ddd(2)
204            enddo
205            xq(2)=r(1,2)*e1(1,i)+r(2,2)*e1(2,i)+r(3,2)*e1(3,i)
206            yq(2)=r(1,2)*e2(1,i)+r(2,2)*e2(2,i)+r(3,2)*e2(3,i)
207!
208            do k=1,3
209               r(k,3)=co(k,kontri(3,j))+vold(k,kontri(3,j))-pmid(k,i)
210            enddo
211            ddd(3)=dsqrt(r(1,3)*r(1,3)+r(2,3)*r(2,3)+r(3,3)*r(3,3))
212            do k=1,3
213               r(k,3)=r(k,3)/ddd(3)
214            enddo
215            xq(3)=r(1,3)*e1(1,i)+r(2,3)*e1(2,i)+r(3,3)*e1(3,i)
216            yq(3)=r(1,3)*e2(1,i)+r(2,3)*e2(2,i)+r(3,3)*e2(3,i)
217!
218c            do l=1,3
219c               do k=1,3
220c                  vertex(k,l)=co(k,kontri(l,j))-pmid(k,i)
221c               enddo
222c               dd=dsqrt(vertex(1,l)**2+vertex(2,l)**2+
223c     &              vertex(3,l)**2)
224c               do k=1,3
225c                  vertex(k,l)=vertex(k,l)/dd
226c               enddo
227c               vertexl(1,l)=vertex(1,l)*e1(1,i)+
228c     &              vertex(2,l)*e1(2,i)+
229c     &              vertex(3,l)*e1(3,i)
230c               vertexl(2,l)=vertex(1,l)*e2(1,i)+
231c     &              vertex(2,l)*e2(2,i)+
232c     &              vertex(3,l)*e2(3,i)
233c            enddo
234!
235!     determining the center of gravity of the projected
236!     triangle
237!
238c            do k=1,2
239c               dirloc(k)=(vertexl(k,1)+vertexl(k,2)+
240c     &              vertexl(k,3))/3.d0
241c            enddo
242            dirloc(1)=(xq(1)+xq(2)+xq(3))/3.d0
243            dirloc(2)=(yq(1)+yq(2)+yq(3))/3.d0
244!
245!     check whether this direction was already covered
246!
247            ix=int((dirloc(1)+1.d0)/dint)+1
248            iy=int((dirloc(2)+1.d0)/dint)+1
249            if(covered(ix,iy).eq.'T') then
250               cycle
251            endif
252!
253!     determine the direction vector in global coordinates
254!
255            do k=1,3
256               dir(k)=(pmid(k,j)-pmid(k,i))/dist(k1)
257            enddo
258!
259!     direction vector in local coordinates of triangle i
260!
261            dirloc(3)=dir(1)*e3(1,i)+dir(2)*e3(2,i)+dir(3)*e3(3,i)
262!
263!     if surfaces are close, numerical integration with
264!     cubtri is performed
265!
266            if(dist(k1).le.factor*dsqrt(area(i))*dirloc(3)) then
267!
268!     vertices of triangle j
269!
270               do k=1,3
271                  do l=1,3
272                     q(l,k)=co(l,kontri(k,j))+vold(l,kontri(k,j))
273                  enddo
274               enddo
275!
276!     formfactor contribution
277!
278               epsrel=0.01d0
279               epsabs=0.d0
280               limev=100
281               nw=239
282               ncalls=0
283!
284!              storing common data into field rdata
285!
286               rdata(1)=q(1,1)
287               rdata(2)=q(1,2)
288               rdata(3)=q(1,3)
289               rdata(4)=q(2,1)
290               rdata(5)=q(2,2)
291               rdata(6)=q(2,3)
292               rdata(7)=q(3,1)
293               rdata(8)=q(3,2)
294               rdata(9)=q(3,3)
295               rdata(10)=e1(1,i)
296               rdata(11)=e1(2,i)
297               rdata(12)=e1(3,i)
298               rdata(13)=e2(1,i)
299               rdata(14)=e2(2,i)
300               rdata(15)=e2(3,i)
301               rdata(16)=e3(1,i)
302               rdata(17)=e3(2,i)
303               rdata(18)=e3(3,i)
304               rdata(19)=porigin(1)
305               rdata(20)=porigin(2)
306               rdata(21)=porigin(3)
307!
308!     max 1000 evaluations for nw=239
309!
310               call cubtri(fform,pl,epsrel,limev,ftij,abserr,ncalls,
311     &              w,nw,idata,rdata,ier)
312               ftij=ftij/2.d0
313            endif
314!
315!     updating the coverage matrix
316!
317c            do k=1,3
318c               r(k,1)=co(k,kontri(1,j))+vold(k,kontri(1,j))-pmid(k,i)
319c            enddo
320c            ddd(1)=dsqrt(r(1,1)*r(1,1)+r(2,1)*r(2,1)+r(3,1)*r(3,1))
321c            do k=1,3
322c               r(k,1)=r(k,1)/ddd(1)
323c            enddo
324c            xq(1)=r(1,1)*e1(1,i)+r(2,1)*e1(2,i)+r(3,1)*e1(3,i)
325c            yq(1)=r(1,1)*e2(1,i)+r(2,1)*e2(2,i)+r(3,1)*e2(3,i)
326c!
327c            do k=1,3
328c               r(k,2)=co(k,kontri(2,j))+vold(k,kontri(2,j))-pmid(k,i)
329c            enddo
330c            ddd(2)=dsqrt(r(1,2)*r(1,2)+r(2,2)*r(2,2)+r(3,2)*r(3,2))
331c            do k=1,3
332c               r(k,2)=r(k,2)/ddd(2)
333c            enddo
334c            xq(2)=r(1,2)*e1(1,i)+r(2,2)*e1(2,i)+r(3,2)*e1(3,i)
335c            yq(2)=r(1,2)*e2(1,i)+r(2,2)*e2(2,i)+r(3,2)*e2(3,i)
336c!
337c            do k=1,3
338c               r(k,3)=co(k,kontri(3,j))+vold(k,kontri(3,j))-pmid(k,i)
339c            enddo
340c            ddd(3)=dsqrt(r(1,3)*r(1,3)+r(2,3)*r(2,3)+r(3,3)*r(3,3))
341c            do k=1,3
342c               r(k,3)=r(k,3)/ddd(3)
343c            enddo
344c            xq(3)=r(1,3)*e1(1,i)+r(2,3)*e1(2,i)+r(3,3)*e1(3,i)
345c            yq(3)=r(1,3)*e2(1,i)+r(2,3)*e2(2,i)+r(3,3)*e2(3,i)
346!
347            if(dabs(xq(2)-xq(1)).lt.1.d-5) xq(2)=xq(1)+1.d-5
348            if(dabs(xq(2)-xq(1)).lt.1.d-5) xq(2)=xq(1)+1.d-5
349!
350!     if the surfaces are far enough away, one-point
351!     integration is used
352!
353            if(dist(k1).gt.factor*dsqrt(area(i))*dirloc(3)) then
354               ftij=0.d0
355               do k=1,3
356                  l=k-1
357                  if(l.lt.1) l=3
358                  xn(1)=r(2,k)*r(3,l)-r(2,l)*r(3,k)
359                  xn(2)=r(3,k)*r(1,l)-r(3,l)*r(1,k)
360                  xn(3)=r(1,k)*r(2,l)-r(1,l)*r(2,k)
361                  xxn=dsqrt(xn(1)**2+xn(2)**2+xn(3)**2)
362!
363!     argument of dacos must have an absolute value
364!     smaller than or equal to 1.d0; due to
365!     round-off the value can slightly exceed one
366!     and has to be cut-off.
367!
368                  argument=
369     &                 r(1,k)*r(1,l)+r(2,k)*r(2,l)+r(3,k)*r(3,l)
370c     &                 (r(1,k)*r(1,l)+r(2,k)*r(2,l)+r(3,k)*r(3,l))/
371c     &                 (ddd(k)*ddd(l))
372                  if(dabs(argument).gt.1.d0) then
373                     if(argument.gt.0.d0) then
374                        argument=1.d0
375                     else
376                        argument=-1.d0
377                     endif
378                  endif
379                  ftij=ftij+
380     &                 (e3(1,i)*xn(1)
381     &                 +e3(2,i)*xn(2)
382     &                 +e3(3,i)*xn(3))/xxn
383     &                 *dacos(argument)
384               enddo
385               ftij=ftij*area(i)/2.d0
386            endif
387!
388            idj=kontri(4,j)
389            i0=0
390            call add_sm_st_as(auview,adview,jqrad,irowrad,
391     &           idi,idj,ftij,i0,i0,nzsradv)
392!
393!     determining maxima and minima
394!
395            xqmin=2.d0
396            xqmax=-2.d0
397            do k=1,3
398               if(xq(k).lt.xqmin) then
399                  xqmin=xq(k)
400                  imin=k
401               endif
402               if(xq(k).gt.xqmax) then
403                  xqmax=xq(k)
404                  imax=k
405               endif
406            enddo
407!
408            if(((imin.eq.1).and.(imax.eq.2)).or.
409     &           ((imin.eq.2).and.(imax.eq.1))) then
410               imid=3
411               xxmid=xq(3)
412            elseif(((imin.eq.2).and.(imax.eq.3)).or.
413     &              ((imin.eq.3).and.(imax.eq.2))) then
414               imid=1
415               xxmid=xq(1)
416            else
417               imid=2
418               xxmid=xq(2)
419            endif
420!
421!     check for equal x-values
422!
423            if(xxmid-xqmin.lt.1.d-5) then
424               xqmin=xqmin-1.d-5
425               xq(imin)=xqmin
426            endif
427            if(xqmax-xxmid.lt.1.d-5) then
428               xqmax=xqmax+1.d-5
429               xq(imax)=xqmax
430            endif
431!
432!     equation of the straight lines connecting the
433!     triangle vertices in the local x-y plane
434!
435            a(1,2)=yq(2)-yq(1)
436            b(1,2)=xq(1)-xq(2)
437            c(1,2)=yq(1)*xq(2)-xq(1)*yq(2)
438!
439            a(2,3)=yq(3)-yq(2)
440            b(2,3)=xq(2)-xq(3)
441            c(2,3)=yq(2)*xq(3)-xq(2)*yq(3)
442!
443            a(3,1)=yq(1)-yq(3)
444            b(3,1)=xq(3)-xq(1)
445            c(3,1)=yq(3)*xq(1)-xq(3)*yq(1)
446!
447            a(2,1)=a(1,2)
448            b(2,1)=b(1,2)
449            c(2,1)=c(1,2)
450            a(3,2)=a(2,3)
451            b(3,2)=b(2,3)
452            c(3,2)=c(2,3)
453            a(1,3)=a(3,1)
454            b(1,3)=b(3,1)
455            c(1,3)=c(3,1)
456!
457            istart=int((xqmin+1.d0+dint/2.d0)/dint)+1
458            iend=int((xxmid+1.d0+dint/2.d0)/dint)
459            do i1=istart,iend
460               x=dint*(i1-0.5d0)-1.d0
461               yqmin=-(a(imin,imid)*x+c(imin,imid))/b(imin,imid)
462               yqmax=-(a(imin,imax)*x+c(imin,imax))/b(imin,imax)
463               if(yqmin.gt.yqmax) then
464                  dummy=yqmin
465                  yqmin=yqmax
466                  yqmax=dummy
467               endif
468               jstart=int((yqmin+1.d0+dint/2.d0)/dint)+1
469               jend=int((yqmax+1.d0+dint/2.d0)/dint)
470               do j1=jstart,jend
471                  covered(i1,j1)='T'
472               enddo
473               ncovered=ncovered+jend-jstart+1
474            enddo
475!
476            istart=int((xxmid+1.d0+dint/2.d0)/dint)+1
477            iend=int((xqmax+1.d0+dint/2.d0)/dint)
478            do i1=istart,iend
479               x=dint*(i1-0.5d0)-1.d0
480               yqmin=-(a(imid,imax)*x+c(imid,imax))/b(imid,imax)
481               yqmax=-(a(imin,imax)*x+c(imin,imax))/b(imin,imax)
482               if(yqmin.gt.yqmax) then
483                  dummy=yqmin
484                  yqmin=yqmax
485                  yqmax=dummy
486               endif
487               jstart=int((yqmin+1.d0+dint/2.d0)/dint)+1
488               jend=int((yqmax+1.d0+dint/2.d0)/dint)
489               do j1=jstart,jend
490                  covered(i1,j1)='T'
491               enddo
492               ncovered=ncovered+jend-jstart+1
493            enddo
494            if(ncovered.eq.ng*ng)exit
495!
496         enddo
497      enddo
498!
499      return
500      end
501!
502!     function to be integrated
503!
504      real*8 function fform(x,y,idata,rdata)
505!
506      implicit none
507!
508      integer k,l,idata(1)
509!
510      real*8 pint(3),ddd(3),xn(3),q(3,3),
511     &   unitvec(3,3),r(3,3),xxn,x,y,porigin(3),rdata(21)
512!
513!     retrieving common data from field rdata
514!
515      q(1,1)=rdata(1)
516      q(1,2)=rdata(2)
517      q(1,3)=rdata(3)
518      q(2,1)=rdata(4)
519      q(2,2)=rdata(5)
520      q(2,3)=rdata(6)
521      q(3,1)=rdata(7)
522      q(3,2)=rdata(8)
523      q(3,3)=rdata(9)
524      unitvec(1,1)=rdata(10)
525      unitvec(2,1)=rdata(11)
526      unitvec(3,1)=rdata(12)
527      unitvec(1,2)=rdata(13)
528      unitvec(2,2)=rdata(14)
529      unitvec(3,2)=rdata(15)
530      unitvec(1,3)=rdata(16)
531      unitvec(2,3)=rdata(17)
532      unitvec(3,3)=rdata(18)
533      porigin(1)=rdata(19)
534      porigin(2)=rdata(20)
535      porigin(3)=rdata(21)
536!
537      do k=1,3
538         pint(k)=porigin(k)+x*unitvec(k,1)+y*unitvec(k,2)
539      enddo
540!
541      do k=1,3
542         r(k,1)=q(k,1)-pint(k)
543      enddo
544      ddd(1)=dsqrt(r(1,1)*r(1,1)+r(2,1)*r(2,1)+r(3,1)*r(3,1))
545!
546      do k=1,3
547         r(k,2)=q(k,2)-pint(k)
548      enddo
549      ddd(2)=dsqrt(r(1,2)*r(1,2)+r(2,2)*r(2,2)+r(3,2)*r(3,2))
550!
551      do k=1,3
552         r(k,3)=q(k,3)-pint(k)
553      enddo
554      ddd(3)=dsqrt(r(1,3)*r(1,3)+r(2,3)*r(2,3)+r(3,3)*r(3,3))
555!
556!     calculating the contribution
557!
558      fform=0.d0
559      do k=1,3
560         l=k-1
561         if(l.lt.1) l=3
562         xn(1)=r(2,k)*r(3,l)-r(2,l)*r(3,k)
563         xn(2)=r(3,k)*r(1,l)-r(3,l)*r(1,k)
564         xn(3)=r(1,k)*r(2,l)-r(1,l)*r(2,k)
565         xxn=dsqrt(xn(1)**2+xn(2)**2+xn(3)**2)
566         fform=fform+
567     &        (unitvec(1,3)*xn(1)
568     &        +unitvec(2,3)*xn(2)
569     &        +unitvec(3,3)*xn(3))/xxn
570     &        *dacos(
571     &        (r(1,k)*r(1,l)+r(2,k)*r(2,l)+r(3,k)*r(3,l))/
572     &        (ddd(k)*ddd(l)))
573      enddo
574!
575      return
576      end
577
578