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 biotsavart(ipkon,kon,lakon,ne,co,qfx,h0,mi,nka,nkb) 20! 21 implicit none 22! 23! calculates the magnetic intensity due to currents in the phi- 24! domain of an electromagnetic calculation 25! 26 character*8 lakon(*) 27! 28 integer ipkon(*),kon(*),ne,i,nka,nkb,mint3d,konl(26), 29 & j,k,indexe,kk,iflag,mi(*),nope,l 30! 31 real*8 co(3,*),qfx(3,mi(1),*),h0(3,*),xl(3,26),r(3),c2, 32 & con(3),pgauss(3),c1,xi,et,ze,xsj,shp(4,20),weight 33! 34 include "gauss.f" 35! 36 c1=1.d0/(16.d0*datan(1.d0)) 37 iflag=2 38! 39 do j=nka,nkb 40! 41 do k=1,3 42 con(k)=co(k,j) 43 enddo 44! 45 do i=1,ne 46 if(ipkon(i).lt.0) cycle 47! 48! currents are supposed to be modeled by shell elements 49! only 50! 51 if(lakon(i)(7:7).ne.'L') cycle 52! 53 if(lakon(i)(4:5).eq.'8R') then 54 mint3d=1 55 nope=8 56 elseif(lakon(i)(4:4).eq.'8') then 57 mint3d=8 58 nope=8 59 elseif(lakon(i)(4:6).eq.'20R') then 60 mint3d=8 61 nope=20 62 elseif(lakon(i)(4:4).eq.'2') then 63 mint3d=27 64 nope=20 65 elseif(lakon(i)(4:5).eq.'15') then 66 mint3d=9 67 nope=15 68 elseif(lakon(i)(4:4).eq.'6') then 69 mint3d=2 70 nope=6 71 endif 72! 73 indexe=ipkon(i) 74! 75 do l=1,nope 76 konl(l)=kon(indexe+l) 77 do k=1,3 78 xl(k,l)=co(k,konl(l)) 79 enddo 80 enddo 81! 82 do kk=1,mint3d 83! 84 if(lakon(i)(4:5).eq.'8R') then 85 xi=gauss3d1(1,kk) 86 et=gauss3d1(2,kk) 87 ze=gauss3d1(3,kk) 88 weight=weight3d1(kk) 89 elseif((lakon(i)(4:4).eq.'8').or. 90 & (lakon(i)(4:6).eq.'20R')) 91 & then 92 xi=gauss3d2(1,kk) 93 et=gauss3d2(2,kk) 94 ze=gauss3d2(3,kk) 95 weight=weight3d2(kk) 96 elseif(lakon(i)(4:4).eq.'2') then 97 xi=gauss3d3(1,kk) 98 et=gauss3d3(2,kk) 99 ze=gauss3d3(3,kk) 100 weight=weight3d3(kk) 101 elseif(lakon(i)(4:5).eq.'15') then 102 xi=gauss3d8(1,kk) 103 et=gauss3d8(2,kk) 104 ze=gauss3d8(3,kk) 105 weight=weight3d8(kk) 106 elseif(lakon(i)(4:4).eq.'6') then 107 xi=gauss3d7(1,kk) 108 et=gauss3d7(2,kk) 109 ze=gauss3d7(3,kk) 110 weight=weight3d7(kk) 111 endif 112! 113! shape functions 114! 115 if(nope.eq.20) then 116 call shape20h(xi,et,ze,xl,xsj,shp,iflag) 117 elseif(nope.eq.8) then 118 call shape8h(xi,et,ze,xl,xsj,shp,iflag) 119 elseif(nope.eq.15) then 120 call shape15w(xi,et,ze,xl,xsj,shp,iflag) 121 elseif(nope.eq.6) then 122 call shape6w(xi,et,ze,xl,xsj,shp,iflag) 123 endif 124! 125! coordinates of the gauss point 126! 127 do k=1,3 128 pgauss(k)=0.d0 129 do l=1,nope 130 pgauss(k)=pgauss(k)+shp(4,l)*xl(k,l) 131 enddo 132 enddo 133! 134! distance from node to gauss point 135! 136 do k=1,3 137 r(k)=con(k)-pgauss(k) 138 enddo 139 140 c2=weight*xsj/((r(1)*r(1)+r(2)*r(2)+r(3)*r(3))**(1.5d0)) 141! 142 h0(1,j)=h0(1,j)+c2* 143 & (qfx(2,kk,i)*r(3)-qfx(3,kk,i)*r(2)) 144 h0(2,j)=h0(2,j)+c2* 145 & (qfx(3,kk,i)*r(1)-qfx(1,kk,i)*r(3)) 146 h0(3,j)=h0(3,j)+c2* 147 & (qfx(1,kk,i)*r(2)-qfx(2,kk,i)*r(1)) 148 enddo 149 enddo 150! 151 do k=1,3 152 h0(k,j)=h0(k,j)*c1 153 enddo 154 enddo 155! 156 return 157 end 158