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      subroutine advecstiff(nope,voldl,ithermal,xl,nelemload,nelemadvec,
20     &  nload,lakon,xload,istep,time,ttime,dtime,sideload,vold,mi,
21     &  xloadold,reltime,nmethod,s,iinc,iponoel,inoel,ielprop,prop,
22     &  ielmat,shcon,nshcon,rhcon,nrhcon,ntmat_,ipkon,kon,cocon,ncocon,
23     &  ipobody,xbody,ibody)
24!
25!     calculates the stiffness of an advective element.
26!     An advective element consists of a face with a forced convection
27!     film condition and a network node
28!
29      implicit none
30!
31      character*8 lakonl,lakon(*)
32      character*20 sideload(*),sideloadl
33!
34      integer nope,i,ithermal(*),j,nelemload(2,*),nelemadvec,nload,id,
35     &  nelem,ig,mint2d,iflag,istep,jltyp,nfield,mi(*),nmethod,k,iinc,
36     &  node,nopes,iponoel(*),inoel(2,*),ielprop(*),ielmat(mi(3),*),
37     &  ipkon(*),
38     &  nshcon(*),nrhcon(*),ntmat_,kon(*),ncocon(2,*),ipobody(2,*),
39     &  ibody(3,*)
40!
41      real*8 tl2(9),voldl(0:mi(2),9),xl(3,9),sinktemp,xi,et,weight,
42     &  xl2(3,8),xsj2(3),shp2(7,9),coords(3),xs2(3,7),dxsj2,areaj,
43     &  temp,xload(2,*),timeend(2),time,ttime,dtime,field,reltime,
44     &  vold(0:mi(2),*),xloadold(2,*),s(60,60),sref,sref2,prop(*),
45     &  shcon(0:3,ntmat_,*),rhcon(0:1,ntmat_,*),cocon(0:6,ntmat_,*),
46     &  xbody(7,*),heatnod,heatfac
47!
48!
49!
50      include "gauss.f"
51!
52      iflag=2
53!
54      timeend(1)=time
55      timeend(2)=ttime+time
56!
57!     number of nodes in the advective face
58!
59      nopes=nope-1
60!
61!     temperature and displacements in the element's nodes
62!
63      do i=1,nope
64         tl2(i)=voldl(0,i)
65      enddo
66      do i=1,nopes
67         if(ithermal(2).eq.2) then
68            do j=1,3
69               xl2(j,i)=xl(j,i)
70            enddo
71         else
72            do j=1,3
73               xl2(j,i)=xl(j,i)+voldl(j,i)
74            enddo
75         endif
76      enddo
77!
78      call nident2(nelemload,nelemadvec,nload,id)
79!
80!     the second entry in nelemload points to the original
81!     film loading
82!
83      id=nelemload(2,id)
84!
85!     number of the original element
86!
87      nelem=nelemload(1,id)
88      lakonl=lakon(nelem)
89      read(sideload(id)(2:2),'(i1)') ig
90!
91!     number of integration points
92!
93      if(lakonl(4:5).eq.'8R') then
94         mint2d=1
95      elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R')) then
96         if((lakonl(7:7).eq.'A').or.(lakonl(7:7).eq.'E')) then
97            mint2d=2
98         else
99            mint2d=4
100         endif
101      elseif(lakonl(4:4).eq.'2') then
102         mint2d=9
103      elseif(lakonl(4:5).eq.'10') then
104         mint2d=3
105      elseif(lakonl(4:4).eq.'4') then
106         mint2d=1
107      elseif(lakonl(4:5).eq.'15') then
108         if(ig.le.2) then
109            mint2d=3
110         else
111            mint2d=4
112         endif
113      elseif(lakonl(4:4).eq.'6') then
114         mint2d=1
115      endif
116!
117      do i=1,mint2d
118!
119!     copying the sink temperature to ensure the same
120!     value in each integration point (sinktemp can be
121!     changed in subroutine film: requirement from the
122!     thermal people)
123!
124         sinktemp=tl2(nope)
125!
126         if((lakonl(4:5).eq.'8R').or.
127     &        ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then
128            xi=gauss2d1(1,i)
129            et=gauss2d1(2,i)
130            weight=weight2d1(i)
131         elseif((lakonl(4:4).eq.'8').or.
132     &           (lakonl(4:6).eq.'20R').or.
133     &           ((lakonl(4:5).eq.'15').and.(nopes.eq.8))) then
134            xi=gauss2d2(1,i)
135            et=gauss2d2(2,i)
136            weight=weight2d2(i)
137         elseif(lakonl(4:4).eq.'2') then
138            xi=gauss2d3(1,i)
139            et=gauss2d3(2,i)
140            weight=weight2d3(i)
141         elseif((lakonl(4:5).eq.'10').or.
142     &           ((lakonl(4:5).eq.'15').and.(nopes.eq.6))) then
143            xi=gauss2d5(1,i)
144            et=gauss2d5(2,i)
145            weight=weight2d5(i)
146         elseif((lakonl(4:4).eq.'4').or.
147     &           ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then
148            xi=gauss2d4(1,i)
149            et=gauss2d4(2,i)
150            weight=weight2d4(i)
151         endif
152!
153         if(nopes.eq.8) then
154            call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
155         elseif(nopes.eq.4) then
156            call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag)
157         elseif(nopes.eq.6) then
158            call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
159         else
160            call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
161         endif
162!
163         dxsj2=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+
164     &        xsj2(3)*xsj2(3))
165         areaj=dxsj2*weight
166!
167         temp=0.d0
168         do j=1,nopes
169            temp=temp+tl2(j)*shp2(4,j)
170         enddo
171!
172!     for nonuniform load: determine the coordinates of the
173!     point (transferred into the user subroutine)
174!
175         if((sideload(id)(3:4).eq.'NU').or.
176     &        (sideload(id)(5:6).eq.'NU')) then
177            do k=1,3
178               coords(k)=0.d0
179               do j=1,nopes
180                  coords(k)=coords(k)+xl2(k,j)*shp2(4,j)
181               enddo
182            enddo
183            read(sideload(id)(2:2),'(i1)') jltyp
184            jltyp=jltyp+10
185            node=nelemload(2,id)
186            sideloadl=sideload(id)
187            sideloadl(1:1)='F'
188            call film(xload(1,id),sinktemp,temp,istep,
189     &           iinc,timeend,nelem,i,coords,jltyp,field,nfield,
190     &           sideloadl,node,areaj,vold,mi,
191     &           ipkon,kon,lakon,iponoel,inoel,ielprop,prop,ielmat,
192     &           shcon,nshcon,rhcon,nrhcon,ntmat_,cocon,ncocon,
193     &           ipobody,xbody,ibody,heatnod,heatfac)
194c            if(nmethod.eq.1) xload(1,id)=xloadold(1,id)+
195c     &           (xload(1,id)-xloadold(1,id))*reltime
196         endif
197!
198         sref=xload(1,id)*areaj
199         shp2(4,nope)=-1.d0
200!
201         do j=1,nope
202            sref2=sref*shp2(4,j)
203            s(nope,j)=s(nope,j)-sref2
204            do k=1,nopes
205               s(k,j)=s(k,j)+sref2*shp2(4,k)
206            enddo
207         enddo
208      enddo
209!
210!     for some axisymmetric and plane strain elements only half the
211!     number of integration points was used
212!
213      if(((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R')).and.
214     &   ((lakonl(7:7).eq.'A').or.(lakonl(7:7).eq.'E'))) then
215         do i=1,nope
216            do j=1,nope
217               s(i,j)=2.d0*s(i,j)
218            enddo
219         enddo
220      endif
221!
222      return
223      end
224
225