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 velinireltoabs(ibody,xbody,cbody,nbody,set, 20 & istartset,iendset,ialset,nset,veold,mi,ipkon,kon,lakon,co, 21 & itreated) 22! 23! assigns the body forces to the elements by use of field ipobody 24! 25 implicit none 26! 27 character*8 lakon(*) 28 character*81 cbody(*),elset,set(*) 29! 30 integer ibody(3,*),i,j,k,l,m,istartset(*),nbody,kon(*),id, 31 & iendset(*),ialset(*),nset,istat,indexe,mi(*),ipkon(*), 32 & node,nope,itreated(*) 33! 34 real*8 xbody(7,*),veold(0:mi(2),*),co(3,*),om,a(3),xn(3),p(3),dd 35! 36 do k=1,nbody 37 elset=cbody(k) 38! 39! check for centrifugal loads 40! 41 if(ibody(1,k).eq.1) then 42 ibody(1,k)=0 43 else 44 cycle 45 endif 46! 47! determine 48! om: rotational speed 49! a: point on axis 50! xn: unit vector on axis 51! 52 om=dsqrt(xbody(1,k)) 53 a(1)=xbody(2,k) 54 a(2)=xbody(3,k) 55 a(3)=xbody(4,k) 56 xn(1)=xbody(5,k) 57 xn(2)=xbody(6,k) 58 xn(3)=xbody(7,k) 59! 60! check whether element number or set name 61! 62 read(elset,'(i21)',iostat=istat) l 63 if(istat.eq.0) then 64! 65! treat element l 66! 67 indexe=ipkon(l) 68! 69! nope is the number of nodes belonging to the element 70! 71 if(lakon(l)(4:5).eq.'20') then 72 nope=20 73 elseif(lakon(l)(4:4).eq.'8') then 74 nope=8 75 elseif(lakon(l)(4:5).eq.'10') then 76 nope=10 77 elseif(lakon(l)(4:4).eq.'4') then 78 nope=4 79 elseif(lakon(l)(4:5).eq.'15') then 80 nope=15 81 elseif(lakon(l)(4:4).eq.'6') then 82 nope=6 83 elseif(lakon(l)(1:2).eq.'ES') then 84! 85! contact elements need not be treated 86! 87 if(lakon(l)(7:7).ne.'C') then 88 nope=ichar(lakon(l)(8:8))-47 89 endif 90 elseif(lakon(l)(1:4).eq.'MASS') then 91 nope=1 92 endif 93! 94 do i=1,nope 95 node=kon(indexe+i) 96 if(itreated(node).eq.1) cycle 97 do m=1,3 98 p(m)=co(m,node)-a(m) 99 enddo 100 dd=p(1)*xn(1)+p(2)*xn(2)+p(3)*xn(3) 101 do m=1,3 102 p(m)=(p(m)-dd*xn(m))*om 103 enddo 104 veold(1,node)=veold(1,node)+ 105 & xn(2)*p(3)-xn(3)*p(2) 106 veold(2,node)=veold(2,node)+ 107 & xn(3)*p(1)-xn(1)*p(3) 108 veold(3,node)=veold(3,node)+ 109 & xn(1)*p(2)-xn(2)*p(1) 110 itreated(node)=1 111 enddo 112 else 113! 114! set name 115! 116c do i=1,nset 117c if(set(i).eq.elset) exit 118c enddo 119 call cident81(set,elset,nset,id) 120 i=nset+1 121 if(id.gt.0) then 122 if(elset.eq.set(id)) then 123 i=id 124 endif 125 endif 126! 127 do j=istartset(i),iendset(i) 128 if(ialset(j).gt.0) then 129 l=ialset(j) 130! 131! treat element l 132! 133 indexe=ipkon(l) 134! 135! nope is the number of nodes belonging to the element 136! 137 if(lakon(l)(4:5).eq.'20') then 138 nope=20 139 elseif(lakon(l)(4:4).eq.'8') then 140 nope=8 141 elseif(lakon(l)(4:5).eq.'10') then 142 nope=10 143 elseif(lakon(l)(4:4).eq.'4') then 144 nope=4 145 elseif(lakon(l)(4:5).eq.'15') then 146 nope=15 147 elseif(lakon(l)(4:4).eq.'6') then 148 nope=6 149 elseif(lakon(l)(1:2).eq.'ES') then 150! 151! contact elements need not be treated 152! 153 if(lakon(l)(7:7).ne.'C') then 154 nope=ichar(lakon(l)(8:8))-47 155 endif 156 elseif(lakon(l)(1:4).eq.'MASS') then 157 nope=1 158 endif 159! 160 do i=1,nope 161 node=kon(indexe+i) 162 if(itreated(node).eq.1) cycle 163 do m=1,3 164 p(m)=co(m,node)-a(m) 165 enddo 166 dd=p(1)*xn(1)+p(2)*xn(2)+p(3)*xn(3) 167 do m=1,3 168 p(m)=(p(m)-dd*xn(m))*om 169 enddo 170 veold(1,node)=veold(1,node)+ 171 & xn(2)*p(3)-xn(3)*p(2) 172 veold(2,node)=veold(2,node)+ 173 & xn(3)*p(1)-xn(1)*p(3) 174 veold(3,node)=veold(3,node)+ 175 & xn(1)*p(2)-xn(2)*p(1) 176 itreated(node)=1 177 enddo 178 else 179 l=ialset(j-2) 180 do 181 l=l-ialset(j) 182 if(l.ge.ialset(j-1)) exit 183! 184! treat element l 185! 186 indexe=ipkon(l) 187! 188! nope is the number of nodes belonging to the element 189! 190 if(lakon(l)(4:5).eq.'20') then 191 nope=20 192 elseif(lakon(l)(4:4).eq.'8') then 193 nope=8 194 elseif(lakon(l)(4:5).eq.'10') then 195 nope=10 196 elseif(lakon(l)(4:4).eq.'4') then 197 nope=4 198 elseif(lakon(l)(4:5).eq.'15') then 199 nope=15 200 elseif(lakon(l)(4:4).eq.'6') then 201 nope=6 202 elseif(lakon(l)(1:2).eq.'ES') then 203! 204! contact elements need not be treated 205! 206 if(lakon(l)(7:7).ne.'C') then 207 nope=ichar(lakon(l)(8:8))-47 208 endif 209 elseif(lakon(l)(1:4).eq.'MASS') then 210 nope=1 211 endif 212! 213 do i=1,nope 214 node=kon(indexe+i) 215 if(itreated(node).eq.1) cycle 216 do m=1,3 217 p(m)=co(m,node)-a(m) 218 enddo 219 dd=p(1)*xn(1)+p(2)*xn(2)+p(3)*xn(3) 220 do m=1,3 221 p(m)=(p(m)-dd*xn(m))*om 222 enddo 223 veold(1,node)=veold(1,node)+ 224 & xn(2)*p(3)-xn(3)*p(2) 225 veold(2,node)=veold(2,node)+ 226 & xn(3)*p(1)-xn(1)*p(3) 227 veold(3,node)=veold(3,node)+ 228 & xn(1)*p(2)-xn(2)*p(1) 229 itreated(node)=1 230 enddo 231 enddo 232 endif 233 enddo 234 endif 235 enddo 236! 237 return 238 end 239 240