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