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 rigidbodys(inpc,textpart,set,istartset,iendset, 20 & ialset,nset,nset_,nalset,nalset_,ipompc,nodempc,coefmpc, 21 & labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,lakon,ipkon,kon,nk,nk_, 22 & nodeboun,ndirboun,ikboun,ilboun,nboun,nboun_,iperturb,ne_, 23 & ctrl,typeboun,istep,istat,n,iline,ipol,inl,ipoinp,inp,co, 24 & ipoinpc,ier) 25! 26! reading the input deck: *RIGID BODY 27! 28 implicit none 29! 30 character*1 typeboun(*),inpc(*) 31 character*8 lakon(*) 32 character*20 labmpc(*) 33 character*81 set(*),elset,noset 34 character*132 textpart(16) 35! 36 integer istartset(*),iendset(*),ialset(*),ipompc(*), 37 & nodempc(3,*), 38 & nset,nset_,nalset,nalset_,nmpc,nmpc_,mpcfree,nk,nk_,ikmpc(*), 39 & ilmpc(*),ipkon(*),kon(*),inoset,ielset,i,node,ielement,id, 40 & indexe,nope,istep,istat,n,irefnode,irotnode,ne_, 41 & j,idof,k,nodeboun(*),ndirboun(*),ikboun(*),ilboun(*), 42 & nboun,nboun_,key,iperturb(*),ipos,iline,ipol,inl,ipoinp(2,*), 43 & inp(3,*),ipoinpc(0:*),jmin,jmax,ier 44! 45 real*8 coefmpc(3,*),ctrl(*),co(3,*) 46! 47 jmin=1 48 jmax=3 49! 50 if(istep.gt.0) then 51 write(*,*) 52 & '*ERROR reading *RIGID BODY: *RIGID BODY should be placed' 53 write(*,*) ' before all step definitions' 54 ier=1 55 return 56 endif 57! 58! the *RIGID BODY option implies a nonlinear geometric 59! calculation 60! 61 if(iperturb(1).eq.1) then 62 write(*,*) '*ERROR reading *RIGID BODY: the *RIGID BODY option' 63 write(*,*) ' cannot be used in a perturbation step' 64 ier=1 65 return 66 endif 67! 68 elset=' 69 &' 70 noset=' 71 &' 72 irefnode=0 73 irotnode=0 74! 75 do i=2,n 76 if(textpart(i)(1:6).eq.'ELSET=') then 77 if(noset(1:1).eq.' ') then 78 elset(1:80)=textpart(i)(7:86) 79 ipos=index(elset,' ') 80 elset(ipos:ipos)='E' 81 else 82 write(*,*) '*ERROR reading *RIGID BODY: either NSET or' 83 write(*,*) ' ELSET can be specified, not both' 84 ier=1 85 return 86 endif 87 elseif(textpart(i)(1:8).eq.'PINNSET=') then 88 if(elset(1:1).eq.' ') then 89 noset(1:80)=textpart(i)(9:88) 90 ipos=index(noset,' ') 91 noset(ipos:ipos)='N' 92 else 93 write(*,*) '*ERROR reading *RIGID BODY: either NSET or' 94 write(*,*) ' ELSET can be specified, not both' 95 ier=1 96 return 97 endif 98 elseif(textpart(i)(1:5).eq.'NSET=') then 99 if(elset(1:1).eq.' ') then 100 noset(1:80)=textpart(i)(6:85) 101 ipos=index(noset,' ') 102 noset(ipos:ipos)='N' 103 else 104 write(*,*) '*ERROR reading *RIGID BODY: either NSET or' 105 write(*,*) ' ELSET can be specified, not both' 106 ier=1 107 return 108 endif 109 elseif(textpart(i)(1:8).eq.'REFNODE=') then 110 read(textpart(i)(9:18),'(i10)',iostat=istat) irefnode 111 if(istat.gt.0) then 112 call inputerror(inpc,ipoinpc,iline, 113 & "*RIGID BODY%",ier) 114 return 115 endif 116 if(irefnode.gt.nk) then 117 write(*,*) '*ERROR reading *RIGID BODY: ref node', 118 & irefnode 119 write(*,*) ' has not been defined' 120 ier=1 121 return 122 endif 123 elseif(textpart(i)(1:8).eq.'ROTNODE=') then 124 read(textpart(i)(9:18),'(i10)',iostat=istat) irotnode 125 if(istat.gt.0) then 126 call inputerror(inpc,ipoinpc,iline, 127 & "*RIGID BODY%",ier) 128 return 129 endif 130 if(irotnode.gt.nk) then 131 write(*,*) '*ERROR reading *RIGID BODY: rot node', 132 & irotnode 133 write(*,*) ' has not been defined' 134 ier=1 135 return 136 endif 137 else 138 write(*,*) 139 & '*WARNING reading *RIGID BODY: parameter not recognized:' 140 write(*,*) ' ', 141 & textpart(i)(1:index(textpart(i),' ')-1) 142 call inputwarning(inpc,ipoinpc,iline, 143 & "*RIGID BODY%") 144 endif 145 enddo 146! 147! check whether a set was defined 148! 149 if((elset(1:1).eq.' ').and. 150 & (noset(1:1).eq.' ')) then 151 write(*,*) '*WARNING reading *RIGID BODY: no set defined' 152 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 153 & ipoinp,inp,ipoinpc) 154 return 155 endif 156! 157 inoset=0 158 ielset=0 159! 160! checking whether the set exists 161! 162 if(noset(1:1).ne.' ') then 163c do i=1,nset 164c if(set(i).eq.noset) then 165c inoset=i 166c exit 167c endif 168c enddo 169 call cident81(set,noset,nset,id) 170 if(id.gt.0) then 171 if(noset.eq.set(id)) then 172 inoset=id 173 endif 174 endif 175 if(inoset.eq.0) then 176 write(*,*) '*WARNING reading *RIGID BODY: node set ',noset 177 write(*,*) ' does not exist' 178 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 179 & ipoinp,inp,ipoinpc) 180 return 181 endif 182 endif 183! 184 if(elset(1:1).ne.' ') then 185c do i=1,nset 186c if(set(i).eq.elset) then 187c ielset=i 188c exit 189c endif 190c enddo 191 call cident81(set,elset,nset,id) 192 if(id.gt.0) then 193 if(elset.eq.set(id)) then 194 ielset=id 195 endif 196 endif 197 if(ielset.eq.0) then 198 write(*,*) '*WARNING reading *RIGID BODY: element set ', 199 & elset 200 write(*,*) ' does not exist' 201 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 202 & ipoinp,inp,ipoinpc) 203 return 204 endif 205 endif 206! 207! check for the existence of irefnode and irotnode; if none were 208! defined, new nodes are generated 209! 210 if(irefnode.eq.0) then 211 nk=nk+1 212 if(nk.gt.nk_) then 213 write(*,*) '*ERROR reading *RIGID BODY: increase nk_' 214 ier=1 215 return 216 endif 217 irefnode=nk 218! 219! default position of the reference node is the origin 220! 221 co(1,nk)=0.d0 222 co(2,nk)=0.d0 223 co(3,nk)=0.d0 224 endif 225! 226 if(irotnode.eq.0) then 227 nk=nk+1 228 if(nk.gt.nk_) then 229 write(*,*) '*ERROR reading *RIGID BODY: increase nk_' 230 ier=1 231 return 232 endif 233 irotnode=nk 234 endif 235! 236! check whether other equations apply to the dependent nodes 237! 238 if(inoset.ne.0) then 239 do i=istartset(inoset),iendset(inoset) 240 node=ialset(i) 241 if(node.gt.nk_) then 242 write(*,*) '*ERROR reading *RIGID BODY: node ',node 243 write(*,*) ' belonging to set ',noset 244 write(*,*) ' has not been defined' 245 ier=1 246 return 247 endif 248 if((node.eq.irefnode).or.(node.eq.irotnode)) cycle 249 do j=1,3 250 idof=8*(node-1)+j 251 call nident(ikmpc,idof,nmpc,id) 252 if(id.gt.0) then 253 if(ikmpc(id).eq.idof) then 254 write(*,*) '*WARNING reading *RIGID BODY: dof ',j 255 write(*,*) ' of node ',node,' belonging' 256 write(*,*) ' to a rigid body is detected' 257 write(*,*) ' on the dependent side of ' 258 write(*,*) ' another equation; no rigid' 259 write(*,*) ' body constrained applied' 260 endif 261 endif 262 enddo 263 enddo 264 endif 265! 266 if(ielset.ne.0) then 267 do i=istartset(ielset),iendset(ielset) 268 ielement=ialset(i) 269 if(ielement.gt.ne_) then 270 write(*,*) '*ERROR reading *RIGID BODY: element ', 271 & ielement 272 write(*,*) ' belonging to set ',elset 273 write(*,*) ' has not been defined' 274 ier=1 275 return 276 endif 277 if(ipkon(ielement).lt.0) cycle 278 indexe=ipkon(ielement) 279 if(lakon(ielement)(4:4).eq.'2') then 280 nope=20 281 elseif(lakon(ielement)(4:4).eq.'8') then 282 nope=8 283 elseif(lakon(ielement)(4:5).eq.'10') then 284 nope=10 285 elseif(lakon(ielement)(4:4).eq.'4') then 286 nope=4 287 elseif(lakon(ielement)(4:5).eq.'15') then 288 nope=15 289 else 290 nope=6 291 endif 292 do k=indexe+1,indexe+nope 293 node=kon(k) 294 if((node.eq.irefnode).or.(node.eq.irotnode)) cycle 295 do j=1,3 296 idof=8*(node-1)+j 297 call nident(ikmpc,idof,nmpc,id) 298 if(id.gt.0) then 299 if(ikmpc(id).eq.idof) then 300 write(*,*)'*WARNING reading *RIGID BODY: dof ', 301 & j,'of node ',node,' belonging to a' 302 write(*,*)' rigid body is detected on th 303 &e dependent side of another' 304 write(*,*)' equation; no rigid body cons 305 &trained applied' 306 endif 307 endif 308 enddo 309 enddo 310 enddo 311 endif 312! 313! generating the equations in basis form 314! 315! node set 316! 317 if(inoset.ne.0) then 318 do i=istartset(inoset),iendset(inoset) 319 node=ialset(i) 320 if(node.gt.0) then 321 if((node.eq.irefnode).or.(node.eq.irotnode)) cycle 322 call rigidmpc(ipompc,nodempc,coefmpc,irefnode,irotnode, 323 & labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,nk,nk_, 324 & nodeboun,ndirboun,ikboun,ilboun,nboun,nboun_,node, 325 & typeboun,co,jmin,jmax) 326 else 327 node=ialset(i-2) 328 do 329 node=node-ialset(i) 330 if(node.ge.ialset(i-1)) exit 331 if((node.eq.irefnode).or.(node.eq.irotnode)) cycle 332 call rigidmpc(ipompc,nodempc,coefmpc,irefnode, 333 & irotnode,labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc, 334 & nk,nk_,nodeboun,ndirboun,ikboun,ilboun,nboun, 335 & nboun_,node,typeboun,co,jmin,jmax) 336 enddo 337 endif 338 enddo 339 endif 340! 341! element set 342! 343 if(ielset.ne.0) then 344 do i=istartset(ielset),iendset(ielset) 345 ielement=ialset(i) 346 if(ielement.gt.0) then 347 if(ipkon(ielement).lt.0) cycle 348 indexe=ipkon(ielement) 349 if(lakon(ielement)(4:4).eq.'2') then 350 nope=20 351 elseif(lakon(ielement)(4:4).eq.'8') then 352 nope=8 353 elseif(lakon(ielement)(4:5).eq.'10') then 354 nope=10 355 elseif(lakon(ielement)(4:4).eq.'4') then 356 nope=4 357 elseif(lakon(ielement)(4:5).eq.'15') then 358 nope=15 359 else 360 nope=6 361 endif 362 do k=indexe+1,indexe+nope 363 node=kon(k) 364 if((node.eq.irefnode).or.(node.eq.irotnode)) cycle 365 call rigidmpc(ipompc,nodempc,coefmpc,irefnode, 366 & irotnode,labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc, 367 & nk,nk_,nodeboun,ndirboun,ikboun,ilboun,nboun, 368 & nboun_,node,typeboun,co,jmin,jmax) 369 enddo 370 else 371 ielement=ialset(i-2) 372 do 373 ielement=ielement-ialset(i) 374 if(ielement.ge.ialset(i-1)) exit 375 if(ipkon(ielement).lt.0) cycle 376 indexe=ipkon(ielement) 377 if(lakon(ielement)(4:4).eq.'2') then 378 nope=20 379 elseif(lakon(ielement)(4:4).eq.'8') then 380 nope=8 381 elseif(lakon(ielement)(4:5).eq.'10') then 382 nope=10 383 elseif(lakon(ielement)(4:4).eq.'4') then 384 nope=4 385 elseif(lakon(ielement)(4:5).eq.'15') then 386 nope=15 387 else 388 nope=6 389 endif 390 do k=indexe+1,indexe+nope 391 node=kon(k) 392 if((node.eq.irefnode).or.(node.eq.irotnode)) cycle 393 call rigidmpc(ipompc,nodempc,coefmpc,irefnode, 394 & irotnode,labmpc,nmpc,nmpc_,mpcfree,ikmpc, 395 & ilmpc,nk,nk_,nodeboun,ndirboun,ikboun,ilboun, 396 & nboun,nboun_,node,typeboun,co,jmin,jmax) 397 enddo 398 enddo 399 endif 400 enddo 401 endif 402! 403 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 404 & ipoinp,inp,ipoinpc) 405! 406 return 407 end 408 409