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