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 radmatrix(ntr,adrad,aurad,bcr,sideload,nelemload,
27     &     xloadact,lakon,vold,ipkon,kon,co,nloadtr,tarea,tenv,physcon,
28     &     erad,adview,auview,ithermal,iinc,iit,fenv,istep,dtime,ttime,
29     &     time,iviewfile,xloadold,reltime,nmethod,mi,iemchange,nam,
30     &     iamload,jqrad,irowrad,nzsrad)
31!
32      implicit none
33!
34      character*8 lakonl,lakon(*)
35      character*20 sideload(*)
36!
37      integer ntr,nelemload(2,*),nope,nopes,mint2d,i,j,k,l,
38     &     node,ifaceq(8,6),ifacet(6,4),iviewfile,mi(*),
39     &     ifacew(8,5),nelem,ig,index,konl(20),iflag,
40     &     ipkon(*),kon(*),nloadtr(*),i1,ithermal(*),iinc,iit,
41     &     istep,jltyp,nfield,nmethod,iemchange,nam,
42     &     iamload(2,*),irowrad(*),jqrad(*),nzsrad
43!
44      real*8 adrad(*),aurad(*),bcr(ntr,1),xloadact(2,*),h(2),
45     &     xl2(3,8),coords(3),dxsj2,temp,xi,et,weight,xsj2(3),
46     &     vold(0:mi(2),*),co(3,*),shp2(7,8),xs2(3,7),
47     &     areamean,tarea(*),tenv(*),erad(*),fenv(*),physcon(*),
48     &     adview(*),auview(*),tl2(8),tvar(2),field,
49     &     dtime,ttime,time,areaj,xloadold(2,*),reltime,
50     &     fform
51!
52      data ifaceq /4,3,2,1,11,10,9,12,
53     &            5,6,7,8,13,14,15,16,
54     &            1,2,6,5,9,18,13,17,
55     &            2,3,7,6,10,19,14,18,
56     &            3,4,8,7,11,20,15,19,
57     &            4,1,5,8,12,17,16,20/
58      data ifacet /1,3,2,7,6,5,
59     &     1,2,4,5,9,8,
60     &     2,3,4,6,10,9,
61     &     1,4,3,8,10,7/
62      data ifacew /1,3,2,9,8,7,0,0,
63     &     4,5,6,10,11,12,0,0,
64     &     1,2,5,4,7,14,10,13,
65     &     2,3,6,5,8,15,11,14,
66     &     4,6,3,1,12,15,9,13/
67      data iflag /2/
68!
69      external fform
70!
71      include "gauss.f"
72!
73      tvar(1)=time
74      tvar(2)=ttime+time
75!
76!     filling acr and bcr
77!
78      do i1=1,ntr
79         i=nloadtr(i1)
80         nelem=nelemload(1,i)
81         lakonl=lakon(nelem)
82!
83!     calculate the mean temperature of the face
84!
85         read(sideload(i)(2:2),'(i1)') ig
86!
87!     number of nodes and integration points in the face
88!
89         if(lakonl(4:4).eq.'2') then
90            nope=20
91            nopes=8
92         elseif(lakonl(4:4).eq.'8') then
93            nope=8
94            nopes=4
95         elseif(lakonl(4:5).eq.'10') then
96            nope=10
97            nopes=6
98         elseif(lakonl(4:4).eq.'4') then
99            nope=4
100            nopes=3
101         elseif(lakonl(4:5).eq.'15') then
102            nope=15
103         else
104            nope=6
105         endif
106!
107         if(lakonl(4:5).eq.'8R') then
108            mint2d=1
109         elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R'))
110     &           then
111            if(lakonl(7:7).eq.'A') then
112               mint2d=2
113            else
114               mint2d=4
115            endif
116         elseif(lakonl(4:4).eq.'2') then
117            mint2d=9
118         elseif(lakonl(4:5).eq.'10') then
119            mint2d=3
120         elseif(lakonl(4:4).eq.'4') then
121            mint2d=1
122         endif
123!
124         if(lakonl(4:4).eq.'6') then
125            mint2d=1
126            if(ig.le.2) then
127               nopes=3
128            else
129               nopes=4
130            endif
131         endif
132         if(lakonl(4:5).eq.'15') then
133            if(ig.le.2) then
134               mint2d=3
135               nopes=6
136            else
137               mint2d=4
138               nopes=8
139            endif
140         endif
141!
142!     connectivity of the element
143!
144         index=ipkon(nelem)
145         if(index.lt.0) then
146            write(*,*) '*ERROR in radmatrix: element ',nelem
147            write(*,*) '       is not defined'
148            call exit(201)
149         endif
150         do k=1,nope
151            konl(k)=kon(index+k)
152         enddo
153!
154!     coordinates of the nodes belonging to the face
155!
156         if((nope.eq.20).or.(nope.eq.8)) then
157            do k=1,nopes
158               tl2(k)=vold(0,konl(ifaceq(k,ig)))
159!
160               do j=1,3
161                  xl2(j,k)=co(j,konl(ifaceq(k,ig)))+
162     &                 vold(j,konl(ifaceq(k,ig)))
163               enddo
164            enddo
165         elseif((nope.eq.10).or.(nope.eq.4)) then
166            do k=1,nopes
167               tl2(k)=vold(0,konl(ifacet(k,ig)))
168               do j=1,3
169                  xl2(j,k)=co(j,konl(ifacet(k,ig)))+
170     &                 vold(j,konl(ifacet(k,ig)))
171               enddo
172            enddo
173         else
174            do k=1,nopes
175               tl2(k)=vold(0,konl(ifacew(k,ig)))
176               do j=1,3
177                  xl2(j,k)=co(j,konl(ifacew(k,ig)))+
178     &                 vold(j,konl(ifacew(k,ig)))
179               enddo
180            enddo
181         endif
182!
183!     integration to obtain the center of gravity and the mean
184!     temperature; radiation coefficient
185!
186         areamean=0.d0
187         tarea(i1)=0.d0
188!
189         read(sideload(i)(2:2),'(i1)') jltyp
190         jltyp=jltyp+10
191         if(sideload(i)(5:6).ne.'NU') then
192            erad(i1)=xloadact(1,i)
193!
194!           if an amplitude was defined for the emissivity it is
195!           assumed that the emissivity changes with the step, so
196!           acr has to be calculated anew in every iteration
197!
198            if(nam.gt.0) then
199               if(iamload(1,i).ne.0) then
200                  iemchange=1
201               endif
202            endif
203         else
204            erad(i1)=0.d0
205         endif
206!
207         do l=1,mint2d
208            if((lakonl(4:5).eq.'8R').or.
209     &           ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then
210               xi=gauss2d1(1,l)
211               et=gauss2d1(2,l)
212               weight=weight2d1(l)
213            elseif((lakonl(4:4).eq.'8').or.
214     &              (lakonl(4:6).eq.'20R').or.
215     &              ((lakonl(4:5).eq.'15').and.(nopes.eq.8))) then
216               xi=gauss2d2(1,l)
217               et=gauss2d2(2,l)
218               weight=weight2d2(l)
219            elseif(lakonl(4:4).eq.'2') then
220               xi=gauss2d3(1,l)
221               et=gauss2d3(2,l)
222               weight=weight2d3(l)
223            elseif((lakonl(4:5).eq.'10').or.
224     &              ((lakonl(4:5).eq.'15').and.(nopes.eq.6))) then
225               xi=gauss2d5(1,l)
226               et=gauss2d5(2,l)
227               weight=weight2d5(l)
228            elseif((lakonl(4:4).eq.'4').or.
229     &              ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then
230               xi=gauss2d4(1,l)
231               et=gauss2d4(2,l)
232               weight=weight2d4(l)
233            endif
234!
235            if(nopes.eq.8) then
236               call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
237            elseif(nopes.eq.4) then
238               call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag)
239            elseif(nopes.eq.6) then
240               call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
241            else
242               call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
243            endif
244!
245            dxsj2=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+
246     &           xsj2(3)*xsj2(3))
247!
248            temp=0.d0
249            do j=1,nopes
250               temp=temp+tl2(j)*shp2(4,j)
251            enddo
252!
253            tarea(i1)=tarea(i1)+temp*dxsj2*weight
254            areamean=areamean+dxsj2*weight
255!
256            if(sideload(i)(5:6).eq.'NU') then
257               areaj=dxsj2*weight
258               do k=1,3
259                  coords(k)=0.d0
260               enddo
261               do j=1,nopes
262                  do k=1,3
263                     coords(k)=coords(k)+xl2(k,j)*shp2(4,j)
264                  enddo
265               enddo
266               call radiate(h(1),tenv(i1),temp,istep,
267     &              iinc,tvar,nelem,l,coords,jltyp,field,nfield,
268     &              sideload(i),node,areaj,vold,mi,iemchange)
269               if(nmethod.eq.1) h(1)=xloadold(1,i)+
270     &              (h(1)-xloadold(1,i))*reltime
271               erad(i1)=erad(i1)+h(1)
272            endif
273!
274         enddo
275         tarea(i1)=tarea(i1)/areamean-physcon(1)
276         if(sideload(i)(5:6).eq.'NU') then
277            erad(i1)=erad(i1)/mint2d
278         endif
279!
280!        updating the right hand side
281!
282         bcr(i1,1)=physcon(2)*(erad(i1)*tarea(i1)**4+
283     &        (1.d0-erad(i1))*fenv(i1)*tenv(i1)**4)
284c         write(*,*) 'radmatrix ',bcr(i1,1)
285!
286      enddo
287!
288!     adrad and aurad is recalculated only if the viewfactors
289!     or the emissivity changed
290!
291      if(((ithermal(1).eq.3).and.(iviewfile.ge.0)).or.
292     &      ((iit.eq.-1).and.(iviewfile.ne.-2)).or.(iemchange.eq.1).or.
293     &      ((iit.eq.0).and.(abs(nmethod).eq.1))) then
294!
295         do i=1,ntr
296            erad(i)=1.d0-erad(i)
297         enddo
298!
299!     diagonal entries
300!
301         do i=1,ntr
302            adrad(i)=1.d0-erad(i)*adview(i)
303         enddo
304!
305!     lower triangular entries
306!
307         do i=1,nzsrad
308            aurad(i)=-erad(irowrad(i))*auview(i)
309         enddo
310!
311!     upper triangular entries
312!
313         do i=1,ntr
314            do j=nzsrad+jqrad(i),nzsrad+jqrad(i+1)-1
315               aurad(j)=-erad(i)*auview(j)
316            enddo
317         enddo
318!
319         do i=1,ntr
320            erad(i)=1.d0-erad(i)
321         enddo
322      endif
323!
324      return
325      end
326