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 jouleheating(ipkon,lakon,kon,co,elcon,nelcon,
20     &     mi,ne,sti,ielmat,nelemload,sideload,xload,nload,nload_,
21     &     iamload,nam,idefload,ncmat_,ntmat_,
22     &     alcon,nalcon,ithermal,vold,t1,nmethod)
23!
24!     determines the effect of Joule heating
25!
26      implicit none
27!
28      character*8 lakon(*)
29      character*20 label,sideload(*)
30!
31      integer ipkon(*),nelem,kon(*),mi(*),nope,indexe,j,k,null,
32     &     mint3d,jj,iflag,ne,nelemload(2,*),iamload(2,*),nload,nload_,
33     &     ielmat(mi(3),*),konl(20),idefload(*),iamplitude,isector,nam,
34     &     one,nelcon(2,*),nalcon(2,*),ithermal(*),i1,ncmat_,ntmat_,
35     &     imat,nmethod
36!
37      real*8 co(3,*),xl(3,20),xi,et,ze,xsj,shp(4,20),weight,xload(2,*),
38     &     sti(6,mi(1),*),alpha(6),heat,elcon(0:ncmat_,ntmat_,*),volume,
39     &     elconloc(21),t1l,alcon(0:6,ntmat_,*),vold(0:mi(2),*),t1(*)
40!
41      include "gauss.f"
42!
43      data iflag /2/
44!
45      null=0
46      one=1
47!
48      do nelem=1,ne
49        if(ipkon(nelem).lt.0) cycle
50!
51!     only elements belonging to the A,V-domain experience
52!     Joule heating
53!
54        if(int(elcon(2,1,ielmat(1,nelem))).ne.2) cycle
55!
56        imat=ielmat(1,nelem)
57        indexe=ipkon(nelem)
58!
59        if(lakon(nelem)(1:5).eq.'C3D8I') then
60          nope=11
61        elseif(lakon(nelem)(4:4).eq.'2') then
62          nope=20
63        elseif(lakon(nelem)(4:4).eq.'8') then
64          nope=8
65        elseif(lakon(nelem)(4:5).eq.'10') then
66          nope=10
67        elseif(lakon(nelem)(4:4).eq.'4') then
68          nope=4
69        elseif(lakon(nelem)(4:5).eq.'15') then
70          nope=15
71        elseif(lakon(nelem)(4:5).eq.'6') then
72          nope=6
73        endif
74!
75        do j=1,nope
76          konl(j)=kon(indexe+j)
77          do k=1,3
78            xl(k,j)=co(k,konl(j))
79          enddo
80        enddo
81!
82        heat=0.d0
83        volume=0.d0
84!
85        if(lakon(nelem)(4:5).eq.'8R') then
86          mint3d=1
87        elseif((lakon(nelem)(4:4).eq.'8').or.
88     &         (lakon(nelem)(4:6).eq.'20R')) then
89          mint3d=8
90        elseif(lakon(nelem)(4:4).eq.'2') then
91          mint3d=27
92        elseif(lakon(nelem)(4:5).eq.'10') then
93          mint3d=4
94        elseif(lakon(nelem)(4:4).eq.'4') then
95          mint3d=1
96        elseif(lakon(nelem)(4:5).eq.'15') then
97          mint3d=9
98        elseif(lakon(nelem)(4:5).eq.'6') then
99          mint3d=2
100        endif
101!
102!     loop over the integration points
103!
104        do jj=1,mint3d
105          if(lakon(nelem)(4:5).eq.'8R') then
106            xi=gauss3d1(1,jj)
107            et=gauss3d1(2,jj)
108            ze=gauss3d1(3,jj)
109            weight=weight3d1(jj)
110          elseif((lakon(nelem)(4:4).eq.'8').or.
111     &           (lakon(nelem)(4:6).eq.'20R'))
112     &           then
113            xi=gauss3d2(1,jj)
114            et=gauss3d2(2,jj)
115            ze=gauss3d2(3,jj)
116            weight=weight3d2(jj)
117          elseif(lakon(nelem)(4:4).eq.'2') then
118            xi=gauss3d3(1,jj)
119            et=gauss3d3(2,jj)
120            ze=gauss3d3(3,jj)
121            weight=weight3d3(jj)
122          elseif(lakon(nelem)(4:5).eq.'10') then
123            xi=gauss3d5(1,jj)
124            et=gauss3d5(2,jj)
125            ze=gauss3d5(3,jj)
126            weight=weight3d5(jj)
127          elseif(lakon(nelem)(4:4).eq.'4') then
128            xi=gauss3d4(1,jj)
129            et=gauss3d4(2,jj)
130            ze=gauss3d4(3,jj)
131            weight=weight3d4(jj)
132          elseif(lakon(nelem)(4:5).eq.'15') then
133            xi=gauss3d8(1,jj)
134            et=gauss3d8(2,jj)
135            ze=gauss3d8(3,jj)
136            weight=weight3d8(jj)
137          else
138            xi=gauss3d7(1,jj)
139            et=gauss3d7(2,jj)
140            ze=gauss3d7(3,jj)
141            weight=weight3d7(jj)
142          endif
143!
144          if(nope.eq.20) then
145            call shape20h(xi,et,ze,xl,xsj,shp,iflag)
146          elseif(nope.eq.8) then
147            call shape8h(xi,et,ze,xl,xsj,shp,iflag)
148          elseif(nope.eq.10) then
149            call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
150          elseif(nope.eq.4) then
151            call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
152          elseif(nope.eq.15) then
153            call shape15w(xi,et,ze,xl,xsj,shp,iflag)
154          else
155            call shape6w(xi,et,ze,xl,xsj,shp,iflag)
156          endif
157!
158!     calculating the temperature
159!
160          t1l=0.d0
161          if(ithermal(1).eq.1) then
162            if(lakon(nelem)(4:5).eq.'8 ') then
163              do i1=1,nope
164                t1l=t1l+t1(konl(i1))/8.d0
165              enddo
166            elseif(lakon(nelem)(4:6).eq.'20 ')then
167              call linscal(t1,konl,nope,jj,t1l,one)
168            elseif(lakon(nelem)(4:6).eq.'10T') then
169              call linscal10(t1,konl,t1l,null,shp)
170            else
171              do i1=1,nope
172                t1l=t1l+shp(4,i1)*t1(konl(i1))
173              enddo
174            endif
175          elseif(ithermal(1).ge.2) then
176            if(lakon(nelem)(4:5).eq.'8 ') then
177              do i1=1,nope
178                t1l=t1l+vold(0,konl(i1))/8.d0
179              enddo
180            elseif(lakon(nelem)(4:6).eq.'20 ')then
181              call linscal(vold,konl,nope,jj,t1l,mi(2))
182            elseif(lakon(nelem)(4:6).eq.'10T') then
183              call linscal10(vold,konl,t1l,mi(2),shp)
184            else
185              do i1=1,nope
186                t1l=t1l+shp(4,i1)*vold(0,konl(i1))
187              enddo
188            endif
189          endif
190!
191!     material data (electric conductivity and
192!     magnetic permeability)
193!
194          call materialdata_em(elcon,nelcon,alcon,nalcon,
195     &         imat,ntmat_,t1l,elconloc,ncmat_,alpha)
196!
197          if(nmethod.ne.2) then
198!
199!           time dependent current:
200!           heat=sigma*E**2
201!
202            heat=heat+weight*xsj*alpha(1)*
203     &           (sti(1,jj,nelem)**2+
204     &            sti(2,jj,nelem)**2+
205     &            sti(3,jj,nelem)**2)
206          else
207!
208!           alernating current:
209!           heat=sigma*(E_real**2+E_imaginary**2)
210!
211            heat=heat+weight*xsj*alpha(1)*
212     &           (sti(1,jj,nelem)**2+
213     &            sti(2,jj,nelem)**2+
214     &            sti(3,jj,nelem)**2+
215     &            sti(1,jj,nelem+ne)**2+
216     &            sti(2,jj,nelem+ne)**2+
217     &            sti(3,jj,nelem+ne)**2)
218          endif
219          volume=volume+weight*xsj
220        enddo
221!
222        heat=heat/volume
223!
224!     adding the Joule heating to the distributed loading
225!
226        label='BF                  '
227        iamplitude=0
228        isector=0
229        call loadadd(nelem,label,heat,nelemload,sideload,
230     &       xload,nload,nload_,iamload,iamplitude,nam,isector,
231     &       idefload)
232      enddo
233!
234      return
235      end
236
237