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 resultsem(co,kon,ipkon,lakon,v,elcon,nelcon,ielmat,
20     &  ntmat_,vini,dtime,matname,mi,ncmat_,nea,neb,sti,alcon,
21     &  nalcon,h0,istartset,iendset,ialset,iactive,fn,eei,iout,nmethod)
22!
23!     calculates the heat flux and the material tangent at the integration
24!     points and the internal concentrated flux at the nodes
25!
26      implicit none
27!
28      character*8 lakon(*),lakonl
29      character*80 amat,matname(*)
30!
31      integer kon(*),konl(26),mi(*),nelcon(2,*),ielmat(mi(3),*),
32     &  ntmat_,ipkon(*),null,three,iflag,mt,i,j,k,m1,kk,i1,m3,indexe,
33     &  nope,imat,mint3d,ncmat_,nea,neb,nalcon(2,*),mm,l,istart,iset,
34     &  isurf,ilength,istartset(*),iendset(*),ialset(*),nopes,m,
35     &  m2,ig,id,iactive(3),konl2(9),ifaceq(8,6),ifacet(6,4),iel,
36     &  ifacew(8,5),mint2d,nfaces,one,iout,nmethod
37!
38      real*8 co(3,*),v(0:mi(2),*),shp(4,26),xl(3,26),vl(0:mi(2),26),
39     &  elcon(0:ncmat_,ntmat_,*),vkl(0:mi(2),3),vini(0:mi(2),*),c1,
40     &  elconloc(21),xi,et,ze,xsj,t1l,dtime,weight,xsj2(3),shp2(7,9),
41     &  vinikl(0:mi(2),3),alpha(6),vl2(0:mi(2),9),xl2(1:3,9),
42     &  h0l(3),al(3),ainil(3),sti(6,mi(1),*),xs2(3,7),phi,
43     &  um,alcon(0:6,ntmat_,*),h0(3,*),vinil(0:mi(2),26),
44     &  fn(0:mi(2),*),eei(6,mi(1),*)
45!
46      include "gauss.f"
47!
48      data ifaceq /4,3,2,1,11,10,9,12,
49     &            5,6,7,8,13,14,15,16,
50     &            1,2,6,5,9,18,13,17,
51     &            2,3,7,6,10,19,14,18,
52     &            3,4,8,7,11,20,15,19,
53     &            4,1,5,8,12,17,16,20/
54      data ifacet /1,3,2,7,6,5,
55     &             1,2,4,5,9,8,
56     &             2,3,4,6,10,9,
57     &             1,4,3,8,10,7/
58      data ifacew /1,3,2,9,8,7,0,0,
59     &             4,5,6,10,11,12,0,0,
60     &             1,2,5,4,7,14,10,13,
61     &             2,3,6,5,8,15,11,14,
62     &             4,6,3,1,12,15,9,13/
63!
64      iflag=3
65      null=0
66      one=1
67      three=3
68!
69      mt=mi(2)+1
70!
71!     calculation of temperatures and thermal flux
72!
73      do i=nea,neb
74!
75         lakonl(1:8)=lakon(i)(1:8)
76!
77         if(ipkon(i).lt.0) cycle
78!
79         imat=ielmat(1,i)
80         amat=matname(imat)
81!
82         indexe=ipkon(i)
83         if(lakonl(4:5).eq.'20') then
84            nope=20
85            nopes=8
86         elseif(lakonl(4:4).eq.'8') then
87            nope=8
88            nopes=4
89         elseif(lakonl(4:5).eq.'10') then
90            nope=10
91            nopes=6
92         elseif(lakonl(4:4).eq.'4') then
93            nope=4
94            nopes=3
95         elseif(lakonl(4:5).eq.'15') then
96            nope=15
97         elseif(lakonl(4:4).eq.'6') then
98            nope=6
99         else
100            cycle
101         endif
102!
103         if(lakonl(4:5).eq.'8R') then
104            mint3d=1
105            mint2d=1
106         elseif((lakonl(4:4).eq.'8').or.
107     &          (lakonl(4:6).eq.'20R')) then
108            mint3d=8
109            mint2d=4
110         elseif(lakonl(4:4).eq.'2') then
111            mint3d=27
112            mint2d=9
113         elseif(lakonl(4:5).eq.'10') then
114            mint3d=4
115            mint2d=3
116         elseif(lakonl(4:4).eq.'4') then
117            mint3d=1
118            mint2d=1
119         elseif(lakonl(4:5).eq.'15') then
120            mint3d=9
121         elseif(lakonl(4:4).eq.'6') then
122            mint3d=2
123         endif
124!
125         do j=1,nope
126            konl(j)=kon(indexe+j)
127            do k=1,3
128               xl(k,j)=co(k,konl(j))
129            enddo
130            do k=1,5
131               vl(k,j)=v(k,konl(j))
132            enddo
133            vinil(4,j)=vini(4,konl(j))
134         enddo
135!
136         do kk=1,mint3d
137            if(lakonl(4:5).eq.'8R') then
138               xi=gauss3d1(1,kk)
139               et=gauss3d1(2,kk)
140               ze=gauss3d1(3,kk)
141               weight=weight3d1(kk)
142            elseif((lakonl(4:4).eq.'8').or.
143     &             (lakonl(4:6).eq.'20R'))
144     &        then
145               xi=gauss3d2(1,kk)
146               et=gauss3d2(2,kk)
147               ze=gauss3d2(3,kk)
148               weight=weight3d2(kk)
149            elseif(lakonl(4:4).eq.'2') then
150               xi=gauss3d3(1,kk)
151               et=gauss3d3(2,kk)
152               ze=gauss3d3(3,kk)
153               weight=weight3d3(kk)
154            elseif(lakonl(4:5).eq.'10') then
155               xi=gauss3d5(1,kk)
156               et=gauss3d5(2,kk)
157               ze=gauss3d5(3,kk)
158               weight=weight3d5(kk)
159            elseif(lakonl(4:4).eq.'4') then
160               xi=gauss3d4(1,kk)
161               et=gauss3d4(2,kk)
162               ze=gauss3d4(3,kk)
163               weight=weight3d4(kk)
164            elseif(lakonl(4:5).eq.'15') then
165               xi=gauss3d8(1,kk)
166               et=gauss3d8(2,kk)
167               ze=gauss3d8(3,kk)
168               weight=weight3d8(kk)
169            elseif(lakonl(4:4).eq.'6') then
170               xi=gauss3d7(1,kk)
171               et=gauss3d7(2,kk)
172               ze=gauss3d7(3,kk)
173               weight=weight3d7(kk)
174            endif
175!
176            if(nope.eq.20) then
177               call shape20h(xi,et,ze,xl,xsj,shp,iflag)
178            elseif(nope.eq.8) then
179               call shape8h(xi,et,ze,xl,xsj,shp,iflag)
180            elseif(nope.eq.10) then
181               call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
182            elseif(nope.eq.4) then
183               call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
184            elseif(nope.eq.15) then
185               call shape15w(xi,et,ze,xl,xsj,shp,iflag)
186            else
187               call shape6w(xi,et,ze,xl,xsj,shp,iflag)
188            endif
189!
190            c1=xsj*weight
191!
192!                 vkl(m2,m3) contains the derivative of the m2-
193!                 component of the displacement with respect to
194!                 direction m3
195!
196            do k=1,5
197               do m3=1,3
198                  vkl(k,m3)=0.d0
199               enddo
200!
201               do m1=1,nope
202                  do m3=1,3
203                     vkl(k,m3)=vkl(k,m3)+shp(m3,m1)*vl(k,m1)
204                  enddo
205               enddo
206            enddo
207!
208            do m3=1,3
209               vinikl(4,m3)=0.d0
210            enddo
211            do m1=1,nope
212               do m3=1,3
213                  vinikl(4,m3)=vinikl(4,m3)+shp(m3,m1)*vinil(4,m1)
214               enddo
215            enddo
216!
217!              calculating the temperature difference in
218!              the integration point
219!
220            t1l=0.d0
221            do j=1,3
222               h0l(j)=0.d0
223               al(j)=0.d0
224               ainil(j)=0.d0
225            enddo
226            if(lakonl(4:5).eq.'8 ') then
227               do i1=1,8
228                  do j=1,3
229                     h0l(j)=h0l(j)+h0(j,konl(i1))/8.d0
230                     al(j)=al(j)+v(j,konl(i1))/8.d0
231                     ainil(j)=ainil(j)+vini(j,konl(i1))/8.d0
232                  enddo
233                  t1l=t1l+v(0,konl(i1))/8.d0
234               enddo
235            elseif(lakonl(4:6).eq.'20 ') then
236               call linscal(v,konl,nope,kk,t1l,mi(2))
237               call linvec(h0,konl,nope,kk,h0l,one,three)
238               call linvec(v,konl,nope,kk,al,null,mi(2))
239               call linvec(vini,konl,nope,kk,ainil,null,mi(2))
240            elseif(lakonl(4:6).eq.'10T') then
241               call linscal10(v,konl,t1l,mi(2),shp)
242               call linvec10(h0,konl,h0l,one,three,shp)
243               call linvec10(v,konl,al,null,mi(2),shp)
244               call linvec10(vini,konl,ainil,null,mi(2),shp)
245            else
246               do i1=1,nope
247                  t1l=t1l+shp(4,i1)*v(0,konl(i1))
248                  do j=1,3
249                     h0l(j)=h0l(j)+shp(4,i1)*h0(j,konl(i1))
250                     al(j)=al(j)+shp(4,i1)*v(j,konl(i1))
251                     ainil(j)=ainil(j)+shp(4,i1)*vini(j,konl(i1))
252                  enddo
253               enddo
254            endif
255!
256!                 material data (permeability)
257!
258            call materialdata_em(elcon,nelcon,alcon,nalcon,
259     &           imat,ntmat_,t1l,elconloc,ncmat_,alpha)
260!
261            um=elconloc(1)
262!
263            if(int(elconloc(2)).eq.1) then
264!
265!              magnetic field B in phi-domain
266!
267               do k=1,3
268                  sti(k+3,kk,i)=um*(h0l(k)-vkl(5,k))
269               enddo
270!
271!     calculating the electromagnetic force K_phiphi
272!
273               do m1=1,nope
274                  do m3=1,3
275                     fn(5,konl(m1))=fn(5,konl(m1))-c1*
276     &                  um*shp(m3,m1)*vkl(5,m3)
277                  enddo
278               enddo
279            else
280!
281!              magnetic field B in A and A-V domain
282!
283               sti(4,kk,i)=vkl(3,2)-vkl(2,3)
284               sti(5,kk,i)=vkl(1,3)-vkl(3,1)
285               sti(6,kk,i)=vkl(2,1)-vkl(1,2)
286!
287!              electric field E in A-V domain
288!
289               if(int(elconloc(2)).eq.2) then
290                 if((nmethod.eq.2).and.(iout.eq.1)) then
291                   do k=1,3
292                     sti(k,kk,i)=-al(k)-vkl(4,k)
293                   enddo
294                 else
295                   do k=1,3
296                     sti(k,kk,i)=(ainil(k)-al(k)+
297     &                    vinikl(4,k)-vkl(4,k))/dtime
298                   enddo
299                 endif
300               endif
301!
302!     calculating the electromagnetic force K_AA
303!
304               do m1=1,nope
305                  do m2=1,3
306                     do m3=1,3
307                        fn(m2,konl(m1))=fn(m2,konl(m1))+c1*
308     &                       (shp(m3,m1)*(vkl(m2,m3)-vkl(m3,m2))+
309     &                       shp(m2,m1)*vkl(m3,m3))/um
310                     enddo
311                  enddo
312               enddo
313!
314            endif
315!
316         enddo
317!
318!        surface integrals
319!
320!        determining the number of faces per element
321!
322         if((lakonl(4:4).eq.'8').or.(lakonl(4:4).eq.'2')) then
323            nfaces=6
324         elseif((lakonl(4:4).eq.'6').or.(lakonl(4:5).eq.'15')) then
325            nfaces=5
326         elseif((lakonl(4:4).eq.'4').or.(lakonl(4:5).eq.'10')) then
327            nfaces=4
328         endif
329!
330         m=int(elconloc(2))
331            if(iactive(m).eq.0) cycle
332            iset=iactive(m)
333            istart=istartset(iset)
334            ilength=iendset(iset)-istart+1
335!
336            isurf=10*i+nfaces
337            call nident(ialset(istart),isurf,ilength,id)
338!
339            do
340               if(id.eq.0) exit
341               isurf=ialset(istart+id-1)
342               iel=int(isurf/10.d0)
343               if(iel.ne.i) exit
344               ig=isurf-10*iel
345!
346!     treatment of wedge faces
347!
348               if(lakonl(4:4).eq.'6') then
349                  mint2d=1
350                  if(ig.le.2) then
351                     nopes=3
352                  else
353                     nopes=4
354                  endif
355               endif
356               if(lakonl(4:5).eq.'15') then
357                  if(ig.le.2) then
358                     mint2d=3
359                     nopes=6
360                  else
361                     mint2d=4
362                     nopes=8
363                  endif
364               endif
365!
366!              face connectivity is stored in konl2
367!
368               if((nope.eq.20).or.(nope.eq.8)) then
369                  do j=1,nopes
370                     konl2(j)=konl(ifaceq(j,ig))
371                  enddo
372               elseif((nope.eq.10).or.(nope.eq.4)) then
373                  do j=1,nopes
374                     konl2(j)=konl(ifacet(j,ig))
375                  enddo
376               else
377                  do j=1,nopes
378                     konl2(j)=konl(ifacew(j,ig))
379                  enddo
380               endif
381!
382!              face coordinates and solution fields
383!
384               do j=1,nopes
385                  do k=1,3
386                     xl2(k,j)=co(k,konl2(j))
387                     vl2(k,j)=v(k,konl2(j))
388                  enddo
389                  vl2(5,j)=v(5,konl2(j))
390               enddo
391!
392               do kk=1,mint2d
393!
394                  if((lakonl(4:5).eq.'8R').or.
395     &                 ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then
396                     xi=gauss2d1(1,kk)
397                     et=gauss2d1(2,kk)
398                     weight=weight2d1(kk)
399                  elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R')
400     &                .or.((lakonl(4:5).eq.'15').and.(nopes.eq.8))) then
401                     xi=gauss2d2(1,kk)
402                     et=gauss2d2(2,kk)
403                     weight=weight2d2(kk)
404                  elseif(lakonl(4:4).eq.'2') then
405                     xi=gauss2d3(1,kk)
406                     et=gauss2d3(2,kk)
407                     weight=weight2d3(kk)
408                  elseif((lakonl(4:5).eq.'10').or.
409     &                    ((lakonl(4:5).eq.'15').and.(nopes.eq.6))) then
410                     xi=gauss2d5(1,kk)
411                     et=gauss2d5(2,kk)
412                     weight=weight2d5(kk)
413                  elseif((lakonl(4:4).eq.'4').or.
414     &                    ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then
415                     xi=gauss2d4(1,kk)
416                     et=gauss2d4(2,kk)
417                     weight=weight2d4(kk)
418                  endif
419!
420                  if(nopes.eq.8) then
421                     call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
422                  elseif(nopes.eq.4) then
423                     call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag)
424                  elseif(nopes.eq.6) then
425                     call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
426                  else
427                     call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
428                  endif
429!
430!                 determining the electromagnetic force
431!
432                  if(m.gt.1) then
433!
434!                    determining the values of phi
435!
436                     phi=0.d0
437                     do m1=1,nopes
438                        phi=phi+shp2(4,m1)*vl2(5,m1)
439                     enddo
440!
441                     do m1=1,nopes
442                        do k=1,3
443                           l=k+1
444                           if(l.gt.3) l=1
445                           mm=l+1
446                           if(mm.gt.3) mm=1
447                           fn(k,konl2(m1))=fn(k,konl2(m1))-
448     &                        (shp2(l,m1)*xsj2(mm)-shp2(mm,m1)*xsj2(l))
449     &                          *phi*weight
450                        enddo
451                     enddo
452                  elseif(m.eq.1) then
453!
454!     determining the derivative of A w.r.t. x, y and z
455!
456                     do k=1,3
457                        do m3=1,3
458                           vkl(k,m3)=0.d0
459                        enddo
460                        do m1=1,nopes
461                           do m3=1,3
462                              vkl(k,m3)=vkl(k,m3)+shp2(m3,m1)*vl2(k,m1)
463                           enddo
464                        enddo
465                     enddo
466!
467                     do m1=1,nopes
468                        do k=1,3
469                           l=k+1
470                           if(l.gt.3) l=1
471                           mm=l+1
472                           if(mm.gt.3) mm=1
473                           fn(5,konl2(m1))=fn(5,konl2(m1))+
474     &                        shp2(4,m1)*(vkl(mm,l)-vkl(l,mm))*xsj2(k)
475     &                        *weight
476                        enddo
477                     enddo
478                  endif
479!
480               enddo
481               id=id-1
482            enddo
483      enddo
484!
485      return
486      end
487