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 e_c3d_em(co,konl,lakonl,s,sm,
20     &  ff,nelem,nmethod,ielmat,
21     &  ntmat_,t1,ithermal,vold,idist,
22     &  matname,mi,mass,rhsi,
23     &  ncmat_,elcon,nelcon,h0,iactive,
24     &  alcon,nalcon,istartset,iendset,ialset)
25!
26!     computation of the element matrix and rhs for the element with
27!     the topology in konl
28!
29!     ff: rhs without temperature and eigenstress contribution
30!
31!     nmethod=0: check for positive Jacobian
32!     nmethod=1: stiffness matrix + right hand side
33!     nmethod=2: stiffness matrix + mass matrix
34!     nmethod=3: static stiffness + buckling stiffness
35!     nmethod=4: stiffness matrix + mass matrix
36!
37      implicit none
38!
39      logical mass,rhsi
40!
41      character*8 lakonl
42      character*80 matname(*),amat
43!
44      integer konl(26),ifaceq(8,6),nelem,nmethod,iactive(3),
45     &  ithermal(*),idist,i,j,k,i1,m,one,ii,jj,id,ipointer,ig,kk,mi(*),
46     &  ielmat(mi(3),*),ntmat_,nope,nopes,imat,mint2d,mint3d,
47     &  ifacet(6,4),nopev,ifacew(8,5),ipointeri,ipointerj,iflag,
48     &  nelcon(2,*),ncmat_,nalcon(2,*),iel,ii1,ilength,istart,iset,
49     &  isurf,jj1,istartset(*),iendset(*),ialset(*),three,nfaces,
50     &  null
51!
52      real*8 co(3,*),xl(3,26),shp(4,26),s(100,100),ff(100),xs2(3,7),
53     &  t1(*),h0(3,*),xl2(3,9),xsj2(3),shp2(7,9),vold(0:mi(2),*),
54     &  xi,et,ze,xsj,sm(100,100),t1l,weight,
55     &  elcon(0:ncmat_,ntmat_,*),elconloc(21),sigma,
56     &  h0l(3),h0l2(3,8),alcon(0:6,ntmat_,*),diag,um,uminv,alpha(6),
57     &  d(3,3)
58!
59      include "gauss.f"
60!
61      data ifaceq /4,3,2,1,11,10,9,12,
62     &            5,6,7,8,13,14,15,16,
63     &            1,2,6,5,9,18,13,17,
64     &            2,3,7,6,10,19,14,18,
65     &            3,4,8,7,11,20,15,19,
66     &            4,1,5,8,12,17,16,20/
67      data ifacet /1,3,2,7,6,5,
68     &             1,2,4,5,9,8,
69     &             2,3,4,6,10,9,
70     &             1,4,3,8,10,7/
71      data ifacew /1,3,2,9,8,7,0,0,
72     &             4,5,6,10,11,12,0,0,
73     &             1,2,5,4,7,14,10,13,
74     &             2,3,6,5,8,15,11,14,
75     &             4,6,3,1,12,15,9,13/
76      data d /1.d0,0.d0,0.d0,0.d0,1.d0,0.d0,0.d0,0.d0,1.d0/
77!
78      null=0
79      one=1
80      three=3
81      iflag=3
82!
83      imat=ielmat(1,nelem)
84      amat=matname(imat)
85!
86      if(lakonl(4:5).eq.'20') then
87         nope=20
88         nopev=8
89         nopes=8
90      elseif(lakonl(4:4).eq.'8') then
91         nope=8
92         nopev=8
93         nopes=4
94      elseif(lakonl(4:5).eq.'10') then
95         nope=10
96         nopev=4
97         nopes=6
98      elseif(lakonl(4:4).eq.'4') then
99         nope=4
100         nopev=4
101         nopes=3
102      elseif(lakonl(4:5).eq.'15') then
103         nope=15
104         nopev=6
105      elseif(lakonl(4:4).eq.'6') then
106         nope=6
107         nopev=6
108      endif
109!
110      if(lakonl(4:5).eq.'8R') then
111         mint2d=1
112         mint3d=1
113      elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R')) then
114         mint2d=4
115         mint3d=8
116      elseif(lakonl(4:4).eq.'2') then
117         mint2d=9
118         mint3d=27
119      elseif(lakonl(4:5).eq.'10') then
120         mint2d=3
121         mint3d=4
122      elseif(lakonl(4:4).eq.'4') then
123         mint2d=1
124         mint3d=1
125      elseif(lakonl(4:5).eq.'15') then
126         mint3d=9
127      elseif(lakonl(4:4).eq.'6') then
128         mint3d=2
129      else
130         mint3d=0
131      endif
132!
133!     computation of the coordinates of the local nodes
134!
135      do i=1,nope
136        do j=1,3
137          xl(j,i)=co(j,konl(i))
138        enddo
139      enddo
140!
141!       initialisation for distributed forces
142!
143      if(rhsi) then
144        if(idist.ne.0) then
145          do i=1,5*nope
146            ff(i)=0.d0
147          enddo
148        endif
149      endif
150!
151!     initialisation of sm
152!
153      if(mass) then
154        do i=1,5*nope
155          do j=i,5*nope
156            sm(i,j)=0.d0
157          enddo
158        enddo
159      endif
160!
161!     initialisation of s
162!
163      do i=1,5*nope
164        do j=i,5*nope
165          s(i,j)=0.d0
166        enddo
167      enddo
168!
169!     computation of the matrix: loop over the Gauss points
170!
171      do kk=1,mint3d
172         if(lakonl(4:5).eq.'8R') then
173            xi=gauss3d1(1,kk)
174            et=gauss3d1(2,kk)
175            ze=gauss3d1(3,kk)
176            weight=weight3d1(kk)
177         elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R'))
178     &           then
179            xi=gauss3d2(1,kk)
180            et=gauss3d2(2,kk)
181            ze=gauss3d2(3,kk)
182            weight=weight3d2(kk)
183         elseif(lakonl(4:4).eq.'2') then
184            xi=gauss3d3(1,kk)
185            et=gauss3d3(2,kk)
186            ze=gauss3d3(3,kk)
187            weight=weight3d3(kk)
188         elseif(lakonl(4:5).eq.'10') then
189            xi=gauss3d5(1,kk)
190            et=gauss3d5(2,kk)
191            ze=gauss3d5(3,kk)
192            weight=weight3d5(kk)
193         elseif(lakonl(4:4).eq.'4') then
194            xi=gauss3d4(1,kk)
195            et=gauss3d4(2,kk)
196            ze=gauss3d4(3,kk)
197            weight=weight3d4(kk)
198         elseif(lakonl(4:5).eq.'15') then
199            xi=gauss3d8(1,kk)
200            et=gauss3d8(2,kk)
201            ze=gauss3d8(3,kk)
202            weight=weight3d8(kk)
203         else
204            xi=gauss3d7(1,kk)
205            et=gauss3d7(2,kk)
206            ze=gauss3d7(3,kk)
207            weight=weight3d7(kk)
208         endif
209!
210!     calculation of the shape functions and their derivatives
211!     in the gauss point
212!
213         if(nope.eq.20) then
214            call shape20h(xi,et,ze,xl,xsj,shp,iflag)
215         elseif(nope.eq.8) then
216            call shape8h(xi,et,ze,xl,xsj,shp,iflag)
217         elseif(nope.eq.10) then
218            call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
219         elseif(nope.eq.4) then
220            call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
221         elseif(nope.eq.15) then
222            call shape15w(xi,et,ze,xl,xsj,shp,iflag)
223         else
224            call shape6w(xi,et,ze,xl,xsj,shp,iflag)
225         endif
226!
227!           check the jacobian determinant
228!
229         if(xsj.lt.1.d-20) then
230            write(*,*) '*ERROR in e_c3d_th: nonpositive jacobian'
231            write(*,*) '       determinant in element',nelem
232            write(*,*)
233            xsj=dabs(xsj)
234            nmethod=0
235         endif
236!
237!        calculating the temperature and the magnetic intensity
238!        in the integration point (one order lower than the
239!        interpolation of the primary unknowns)
240!
241         do j=1,3
242            h0l(j)=0.d0
243         enddo
244         t1l=0.d0
245!
246!        calculating the magnetic intensity
247!
248         if(lakonl(4:5).eq.'8 ') then
249            do i1=1,nope
250               do j=1,3
251                  h0l(j)=h0l(j)+h0(j,konl(i1))/8.d0
252               enddo
253            enddo
254         elseif(lakonl(4:6).eq.'20 ') then
255            call linvec(h0,konl,nope,kk,h0l,one,three)
256         elseif(lakonl(4:6).eq.'10T') then
257            call linvec10(h0,konl,h0l,one,three,shp)
258         else
259            do i1=1,nope
260               do j=1,3
261                  h0l(j)=h0l(j)+shp(4,i1)*h0(j,konl(i1))
262               enddo
263            enddo
264         endif
265!
266!        calculating the temperature
267!
268         if(ithermal(1).eq.1) then
269            if(lakonl(4:5).eq.'8 ') then
270               do i1=1,nope
271                  t1l=t1l+t1(konl(i1))/8.d0
272               enddo
273            elseif(lakonl(4:6).eq.'20 ')then
274               call linscal(t1,konl,nope,kk,t1l,one)
275            elseif(lakonl(4:6).eq.'10T') then
276               call linscal10(t1,konl,t1l,null,shp)
277            else
278               do i1=1,nope
279                  t1l=t1l+shp(4,i1)*t1(konl(i1))
280               enddo
281            endif
282         elseif(ithermal(1).ge.2) then
283            if(lakonl(4:5).eq.'8 ') then
284               do i1=1,nope
285                  t1l=t1l+vold(0,konl(i1))/8.d0
286               enddo
287            elseif(lakonl(4:6).eq.'20 ')then
288               call linscal(vold,konl,nope,kk,t1l,mi(2))
289            elseif(lakonl(4:6).eq.'10T') then
290               call linscal10(vold,konl,t1l,mi(2),shp)
291            else
292               do i1=1,nope
293                  t1l=t1l+shp(4,i1)*vold(0,konl(i1))
294               enddo
295            endif
296         endif
297!
298!        material data (electric conductivity and
299!        magnetic permeability)
300!
301         call materialdata_em(elcon,nelcon,alcon,nalcon,
302     &        imat,ntmat_,t1l,elconloc,ncmat_,alpha)
303!
304         weight=weight*xsj
305!
306         uminv=weight/elconloc(1)
307         um=weight*elconloc(1)
308         sigma=weight*alpha(1)
309!
310         jj1=0
311         do jj=1,nope
312            ii1=0
313            do ii=1,jj
314!
315               if((int(elconloc(2)).eq.2).or.(int(elconloc(2)).eq.3))
316     &             then
317!
318!                 K_AA matrix
319!
320                  diag=shp(1,ii)*shp(1,jj)+shp(2,ii)*shp(2,jj)+
321     &                 shp(3,ii)*shp(3,jj)
322                  do k=1,3
323                     do m=1,3
324                        s(ii1+k,jj1+m)=s(ii1+k,jj1+m)+
325     &                    (diag*d(k,m)-shp(m,ii)*shp(k,jj)
326     &                                +shp(k,ii)*shp(m,jj))*uminv
327                     enddo
328                  enddo
329!
330!                 M_AA matrix
331!
332                  if(mass) then
333                     do k=1,3
334                        sm(ii1+k,jj1+k)=sm(ii1+k,jj1+k)+
335     &                       sigma*shp(4,ii)*shp(4,jj)
336                     enddo
337                  endif
338                  if((int(elconloc(2)).eq.2).and.mass) then
339!
340!                    M_AV and M_VA matrix
341!
342                     do k=1,3
343                        sm(ii1+k,jj1+4)=sm(ii1+k,jj1+4)+
344     &                       sigma*shp(4,ii)*shp(k,jj)
345                        sm(ii1+4,jj1+k)=sm(ii1+4,jj1+k)+
346     &                       sigma*shp(k,ii)*shp(4,jj)
347                     enddo
348!
349!                    M_VV matrix
350!
351                     sm(ii1+4,jj1+4)=sm(ii1+4,jj1+4)+
352     &                    sigma*(shp(1,ii)*shp(1,jj)+
353     &                    shp(2,ii)*shp(2,jj)+
354     &                    shp(3,ii)*shp(3,jj))
355                  endif
356               else if(int(elconloc(2)).eq.1) then
357!
358!                 K_phiphi matrix
359!
360                  s(ii1+5,jj1+5)=s(ii1+5,jj1+5)-
361     &                     um*(shp(1,ii)*shp(1,jj)+
362     &                    shp(2,ii)*shp(2,jj)+
363     &                    shp(3,ii)*shp(3,jj))
364               endif
365!
366               ii1=ii1+5
367            enddo
368!
369!           F_phi matrix
370!
371            if((rhsi).and.(int(elconloc(2)).eq.1)) then
372               ff(jj1+5)=ff(jj1+5)-um*(shp(1,jj)*h0l(1)+
373     &                                 shp(2,jj)*h0l(2)+
374     &                                 shp(3,jj)*h0l(3))
375            endif
376!
377            jj1=jj1+5
378         enddo
379      enddo
380!
381!     surface integrals
382!
383!     determining the number of faces per element
384!
385      if((lakonl(4:4).eq.'8').or.(lakonl(4:4).eq.'2')) then
386         nfaces=6
387      elseif((lakonl(4:4).eq.'6').or.(lakonl(4:5).eq.'15')) then
388         nfaces=5
389      elseif((lakonl(4:4).eq.'4').or.(lakonl(4:5).eq.'10')) then
390         nfaces=4
391      endif
392!
393      m=int(elconloc(2))
394      if(iactive(m).eq.0) return
395      iset=iactive(m)
396      istart=istartset(iset)
397      ilength=iendset(iset)-istart+1
398!
399      isurf=10*nelem+nfaces
400      call nident(ialset(istart),isurf,ilength,id)
401!
402      do
403         if(id.eq.0) exit
404         isurf=ialset(istart+id-1)
405         iel=int(isurf/10.d0)
406         if(iel.ne.nelem) exit
407         ig=isurf-10*iel
408!
409!     treatment of wedge faces
410!
411         if(lakonl(4:4).eq.'6') then
412            mint2d=1
413            if(ig.le.2) then
414               nopes=3
415            else
416               nopes=4
417            endif
418         endif
419         if(lakonl(4:5).eq.'15') then
420            if(ig.le.2) then
421               mint2d=3
422               nopes=6
423            else
424               mint2d=4
425               nopes=8
426            endif
427         endif
428!
429         if((nope.eq.20).or.(nope.eq.8)) then
430            do i=1,nopes
431               do j=1,3
432                  h0l2(j,i)=h0(j,konl(ifaceq(i,ig)))
433                  xl2(j,i)=co(j,konl(ifaceq(i,ig)))
434               enddo
435            enddo
436         elseif((nope.eq.10).or.(nope.eq.4)) then
437            do i=1,nopes
438               do j=1,3
439                  h0l2(j,i)=h0(j,konl(ifacet(i,ig)))
440                  xl2(j,i)=co(j,konl(ifacet(i,ig)))
441               enddo
442            enddo
443         else
444            do i=1,nopes
445               do j=1,3
446                  h0l2(j,i)=h0(j,konl(ifacew(i,ig)))
447                  xl2(j,i)=co(j,konl(ifacew(i,ig)))
448               enddo
449            enddo
450         endif
451!
452         do i=1,mint2d
453!
454            if((lakonl(4:5).eq.'8R').or.
455     &           ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then
456               xi=gauss2d1(1,i)
457               et=gauss2d1(2,i)
458               weight=weight2d1(i)
459            elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R').or.
460     &              ((lakonl(4:5).eq.'15').and.(nopes.eq.8))) then
461               xi=gauss2d2(1,i)
462               et=gauss2d2(2,i)
463               weight=weight2d2(i)
464            elseif(lakonl(4:4).eq.'2') then
465               xi=gauss2d3(1,i)
466               et=gauss2d3(2,i)
467               weight=weight2d3(i)
468            elseif((lakonl(4:5).eq.'10').or.
469     &              ((lakonl(4:5).eq.'15').and.(nopes.eq.6))) then
470               xi=gauss2d5(1,i)
471               et=gauss2d5(2,i)
472               weight=weight2d5(i)
473            elseif((lakonl(4:4).eq.'4').or.
474     &              ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then
475               xi=gauss2d4(1,i)
476               et=gauss2d4(2,i)
477               weight=weight2d4(i)
478            endif
479!
480            if(nopes.eq.8) then
481               call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
482            elseif(nopes.eq.4) then
483               call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag)
484            elseif(nopes.eq.6) then
485               call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
486            else
487               call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
488            endif
489!
490            do k=1,3
491               h0l(k)=0.d0
492               do j=1,nopes
493                  h0l(k)=h0l(k)+h0l2(k,j)*shp2(4,j)
494               enddo
495            enddo
496!
497            if((rhsi).and.(m.gt.1)) then
498               do k=1,nopes
499                  if((nope.eq.20).or.(nope.eq.8)) then
500                     ipointer=5*(ifaceq(k,ig)-1)
501                  elseif((nope.eq.10).or.(nope.eq.4)) then
502                     ipointer=5*(ifacet(k,ig)-1)
503                  else
504                     ipointer=5*(ifacew(k,ig)-1)
505                  endif
506!
507!                    F_A vector
508!
509                  ff(ipointer+1)=ff(ipointer+1)+
510     &                 shp2(4,k)*(h0l(2)*xsj2(3)-h0l(3)*xsj2(2))*weight
511                  ff(ipointer+2)=ff(ipointer+2)+
512     &                 shp2(4,k)*(h0l(3)*xsj2(1)-h0l(1)*xsj2(3))*weight
513                  ff(ipointer+3)=ff(ipointer+3)+
514     &                 shp2(4,k)*(h0l(1)*xsj2(2)-h0l(2)*xsj2(1))*weight
515!
516               enddo
517            endif
518!
519            do ii=1,nopes
520               if((nope.eq.20).or.(nope.eq.8)) then
521                  ipointeri=5*(ifaceq(ii,ig)-1)
522               elseif((nope.eq.10).or.(nope.eq.4)) then
523                  ipointeri=5*(ifacet(ii,ig)-1)
524               else
525                  ipointeri=5*(ifacew(ii,ig)-1)
526               endif
527               do jj=1,nopes
528                  if((nope.eq.20).or.(nope.eq.8)) then
529                     ipointerj=5*(ifaceq(jj,ig)-1)
530                  elseif((nope.eq.10).or.(nope.eq.4)) then
531                     ipointerj=5*(ifacet(jj,ig)-1)
532                  else
533                     ipointerj=5*(ifacew(jj,ig)-1)
534                  endif
535!
536                  if(m.gt.1) then
537!
538!                       K_Aphi matrix
539!
540                     if((ipointeri+1).gt.ipointerj+5) cycle
541                     s(ipointeri+1,ipointerj+5)=
542     &                    s(ipointeri+1,ipointerj+5)
543     &                    +shp2(4,jj)*(shp2(3,ii)*xsj2(2)
544     &                    -shp2(2,ii)*xsj2(3))
545     &                    *weight
546!
547                     if((ipointeri+2).gt.ipointerj+5) cycle
548                     s(ipointeri+2,ipointerj+5)=
549     &                    s(ipointeri+2,ipointerj+5)
550     &                    +shp2(4,jj)*(shp2(1,ii)*xsj2(3)
551     &                    -shp2(3,ii)*xsj2(1))
552     &                    *weight
553!
554                     if((ipointeri+3).gt.ipointerj+5) cycle
555                     s(ipointeri+3,ipointerj+5)=
556     &                    s(ipointeri+3,ipointerj+5)
557     &                    +shp2(4,jj)*(shp2(2,ii)*xsj2(1)
558     &                    -shp2(1,ii)*xsj2(2))
559     &                    *weight
560!
561!                       K_phiA matrix
562!
563                     if(ipointeri+5.gt.(ipointerj+1)) cycle
564                     s(ipointeri+5,ipointerj+1)=
565     &                    s(ipointeri+5,ipointerj+1)
566     &                    +shp2(4,ii)*(shp2(3,jj)*xsj2(2)
567     &                    -shp2(2,jj)*xsj2(3))
568     &                    *weight
569!
570                     if(ipointeri+5.gt.(ipointerj+2)) cycle
571                     s(ipointeri+5,ipointerj+2)=
572     &                    s(ipointeri+5,ipointerj+2)
573     &                    +shp2(4,ii)*(shp2(1,jj)*xsj2(3)
574     &                    -shp2(3,jj)*xsj2(1))
575     &                    *weight
576!
577                     if(ipointeri+5.gt.(ipointerj+3)) cycle
578                     s(ipointeri+5,ipointerj+3)=
579     &                    s(ipointeri+5,ipointerj+3)
580     &                    +shp2(4,ii)*(shp2(2,jj)*xsj2(1)
581     &                    -shp2(1,jj)*xsj2(2))
582     &                    *weight
583                  endif
584               enddo
585            enddo
586!
587         enddo
588!
589         id=id-1
590      enddo
591!
592      return
593      end
594
595
596
597