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