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