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 geomview(vold,co,pmid,e1,e2,e3,kontri,area, 27 & cs,mcs,inocs,ntrit,nk,mi,sidemean) 28! 29 implicit none 30! 31! change following line if nlabel is increased 32! 33 character*87 label(55) 34! 35 integer i,j,l,mi(*),kontri(4,*),i1,mcs,inocs(*),i2,i3, 36 & ntrit,jj,is,m,nkt,icntrl,imag,nk,nlabel 37! 38 real*8 vold(0:mi(2),*),co(3,*), 39 & pmid(3,*),e3(4,*),e1(3,*),e2(3,*),p1(3),p2(3),p3(3), 40 & cs(17,*),area(*),dd,p21(3),p32(3), 41 & fn,stn,qfn,een,t(3),sidemean,emn 42! 43 write(*,*) 'Calculating the viewfactors' 44 write(*,*) 45! 46! change following line if nlabel is increased and the dimension 47! of field label above! 48! 49 nlabel=55 50! 51! updating the displacements for cyclic symmetric structures 52! 53 if(mcs.gt.0) then 54 nkt=0 55 do i=1,mcs 56 if(int(cs(1,i)).gt.nkt) nkt=int(cs(1,i)) 57 enddo 58 nkt=nk*nkt 59 do i=1,nlabel 60 do l=1,87 61 label(i)(l:l)=' ' 62 enddo 63 enddo 64 label(1)(1:1)='U' 65 imag=0 66 icntrl=2 67 call rectcyl(co,vold,fn,stn,qfn,een,cs,nk,icntrl,t, 68 & label,imag,mi,emn) 69 70 do jj=0,mcs-1 71 is=int(cs(1,jj+1)) 72c write(*,*) 'geomview ',cs(1,jj+1) 73c is=cs(1,jj+1) 74! 75 do i=1,is-1 76 do l=1,nk 77 if(inocs(l).ne.jj) cycle 78 do m=1,mi(2) 79 vold(m,l+nk*i)=vold(m,l) 80 enddo 81 enddo 82 enddo 83 enddo 84 icntrl=-2 85 call rectcyl(co,vold,fn,stn,qfn,een,cs,nkt,icntrl,t, 86 & label,imag,mi,emn) 87 endif 88! 89! calculating the momentaneous center of the triangles, 90! area of the triangles and normal to the triangles 91! 92 sidemean=0.d0 93 do i=1,ntrit 94 i1=kontri(1,i) 95 if(i1.eq.0) cycle 96 i2=kontri(2,i) 97 i3=kontri(3,i) 98 do j=1,3 99 p1(j)=co(j,i1)+vold(j,i1) 100 p2(j)=co(j,i2)+vold(j,i2) 101 p3(j)=co(j,i3)+vold(j,i3) 102 pmid(j,i)=(p1(j)+p2(j)+p3(j))/3.d0 103 p21(j)=p2(j)-p1(j) 104 p32(j)=p3(j)-p2(j) 105 enddo 106! 107! normal to the triangle 108! 109 e3(1,i)=p21(2)*p32(3)-p32(2)*p21(3) 110 e3(2,i)=p21(3)*p32(1)-p32(3)*p21(1) 111 e3(3,i)=p21(1)*p32(2)-p32(1)*p21(2) 112! 113 dd=dsqrt(e3(1,i)*e3(1,i)+e3(2,i)*e3(2,i)+ 114 & e3(3,i)*e3(3,i)) 115! 116! check for degenerated triangles 117! 118 if(dd.lt.1.d-20) then 119 area(i)=0.d0 120 cycle 121 endif 122! 123 do j=1,3 124 e3(j,i)=e3(j,i)/dd 125 enddo 126! 127! area of the triangle 128! 129 area(i)=dd/2.d0 130! 131! unit vector parallel to side 1-2 132! 133 dd=dsqrt(p21(1)*p21(1)+p21(2)*p21(2)+p21(3)*p21(3)) 134 sidemean=sidemean+dd 135 do j=1,3 136 e1(j,i)=p21(j)/dd 137 enddo 138! 139! unit vector orthogonal to e1 and e3 140! 141 e2(1,i)=e3(2,i)*e1(3,i)-e3(3,i)*e1(2,i) 142 e2(2,i)=e3(3,i)*e1(1,i)-e3(1,i)*e1(3,i) 143 e2(3,i)=e3(1,i)*e1(2,i)-e3(2,i)*e1(1,i) 144! 145! the fourth component in e3 is the constant term in the 146! equation of the triangle plane in the form 147! e3(1)*x+e3(2)*y+e3(3)*z+e3(4)=0 148! 149 e3(4,i)=-(e3(1,i)*p1(1)+e3(2,i)*p1(2) 150 & +e3(3,i)*p1(3)) 151 enddo 152 sidemean=sidemean/ntrit 153! 154 return 155 end 156 157 158