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 selectcyclicsymmetrymodess(inpc,textpart,cs,ics, 20 & tieset,istartset, 21 & iendset,ialset,ipompc,nodempc,coefmpc,nmpc,nmpc_,ikmpc,ilmpc, 22 & mpcfree,mcs,set,nset,labmpc,istep,istat,n,iline,ipol,inl, 23 & ipoinp,inp,nmethod,key,ipoinpc,ier) 24! 25! reading the input deck: *SELECT CYCLIC SYMMETRY MODES 26! 27 implicit none 28! 29 character*1 inpc(*) 30 character*20 labmpc(*) 31 character*81 set(*),leftset,tieset(3,*) 32 character*132 textpart(16) 33! 34 integer istep,istat,n,key,i,ns(5),ics(*),istartset(*),ier, 35 & iendset(*),ialset(*),id,ipompc(*),nodempc(3,*),nmpc,nmpc_, 36 & ikmpc(*),ilmpc(*),mpcfree,i1(2),i2(2),i3,i4,i5,j,k, 37 & mpcfreeold,idof,node,ileft,nset,irepeat,ipoinpc(0:*), 38 & mpc,iline,ipol,inl,ipoinp(2,*),inp(3,*),mcs,lprev,ij,nmethod 39! 40 real*8 coefmpc(*),csab(7),x1(2),x2(2),x3,x4,x5,dd,xn,yn,zn, 41 & cs(17,*) 42! 43! irepeat indicates whether the step was preceded by another 44! cyclic symmetry step (irepeat=1) or not (irepeat=0) 45! 46 data irepeat /0/ 47 save irepeat 48! 49 if(istep.eq.0) then 50 write(*,*)'*ERROR reading *SELECT CYCLIC SYMMETRY MODES:' 51 write(*,*)' *SELECT CYCLIC SYMMETRY MODES' 52 write(*,*)' should be placed within a step definition' 53 ier=1 54 return 55 endif 56! 57! check whether in case of cyclic symmetry the frequency procedure 58! is chosen 59! 60 if((nmethod.ne.2).and.(nmethod.ne.13)) then 61 write(*,*) '*ERROR reading *SELECT CYCLIC SYMMETRY MODES:' 62 write(*,*) ' the only valid procedures' 63 write(*,*) ' for cyclic symmetry calculations' 64 write(*,*) ' with nodal diameters are *FREQUENCY' 65 write(*,*) ' and *GREEN' 66 ier=1 67 return 68 endif 69! 70 ns(2)=0 71 ns(3)=0 72! 73 do i=2,n 74 if(textpart(i)(1:5).eq.'NMIN=') then 75 read(textpart(i)(6:15),'(i10)',iostat=istat) ns(2) 76 if(istat.gt.0) then 77 call inputerror(inpc,ipoinpc,iline, 78 & "*SELECT CYCLIC SYMMETRY MODES%",ier) 79 return 80 endif 81 elseif(textpart(i)(1:5).eq.'NMAX=') then 82 read(textpart(i)(6:15),'(i10)',iostat=istat) ns(3) 83 if(istat.gt.0) then 84 call inputerror(inpc,ipoinpc,iline, 85 & "*SELECT CYCLIC SYMMETRY MODES%",ier) 86 return 87 endif 88 else 89 write(*,*) '*WARNING reading *SELECT CYCLIC SYMMETRY MODES:' 90 write(*,*) ' parameter not recognized:' 91 write(*,*) ' ', 92 & textpart(i)(1:index(textpart(i),' ')-1) 93 call inputwarning(inpc,ipoinpc,iline, 94 &"*SELECT CYCLIC SYMMETRY MODES%") 95 endif 96 enddo 97! 98! check the input 99! 100 if(ns(2).lt.0) then 101 ns(2)=0 102 write(*,*) '*WARNING reading *SELECT CYCLIC SYMMETRY MODES:' 103 write(*,*) ' minimum nodal' 104 write(*,*) ' diameter must be nonnegative' 105 endif 106 if(ns(3).lt.ns(2)) then 107 write(*,*) '*ERROR reading *SELECT CYCLIC SYMMETRY MODES:' 108 write(*,*) ' maximum nodal' 109 write(*,*) ' diameter should not exceed minimal one' 110 ier=1 111 return 112 endif 113! 114! loop over all cyclic symmetry parts 115! 116 do ij=1,mcs 117 ns(1)=int(cs(1,ij)) 118 ns(4)=int(cs(4,ij)) 119 leftset=tieset(2,int(cs(17,ij))) 120 lprev=int(cs(14,ij)) 121 do i=1,7 122 csab(i)=cs(5+i,ij) 123 enddo 124! 125! check whether cyclic symmetry axis is part of the structure 126! 127c do i=1,nset 128c if(set(i).eq.leftset) exit 129c enddo 130c ileft=i 131 call cident81(set,leftset,nset,id) 132 ileft=nset+1 133 if(id.gt.0) then 134 if(leftset.eq.set(id)) then 135 ileft=id 136 endif 137 endif 138! 139! if this step was preceded by a cyclic symmetry step: 140! check for MPC's for nodes on the cyclic symmetry axis 141! and delete them 142! 143 if(irepeat.eq.1) then 144 do i=1,ns(4) 145 node=ics(lprev+i) 146 if(node.lt.0) then 147 node=-node 148 do k=1,3 149 idof=8*(node-1)+k 150 call nident(ikmpc,idof,nmpc,id) 151 if(id.gt.0) then 152 if(ikmpc(id).eq.idof) then 153c write(*,*) 'removing MPC',node,k 154 mpc=ilmpc(id) 155 call mpcrem(mpc,mpcfree,nodempc,nmpc,ikmpc, 156 & ilmpc,labmpc,coefmpc,ipompc) 157 endif 158 endif 159 enddo 160 endif 161 enddo 162 endif 163! 164 do i=1,ns(4) 165 node=ics(lprev+i) 166 if(node.lt.0) then 167 node=-node 168 if(ns(2).ne.ns(3)) then 169 if((ns(2).eq.0).or.(ns(2).eq.1)) then 170 write(*,*) '*ERROR: axis of cyclic symmetry' 171 write(*,*) ' is part of the structure;' 172 write(*,*) ' nodal diameters 0, 1, and' 173 write(*,*) ' those above must be each in' 174 write(*,*) ' separate steps.' 175 ier=1 176 return 177 endif 178 endif 179! 180! specifying special MPC's for nodes on the axis 181! 182! normal along the axis 183! 184 xn=csab(4)-csab(1) 185 yn=csab(5)-csab(2) 186 zn=csab(6)-csab(3) 187 dd=dsqrt(xn*xn+yn*yn+zn*zn) 188 xn=xn/dd 189 yn=yn/dd 190 zn=zn/dd 191! 192! nodal diameter 0 193! 194 if(ns(2).eq.0) then 195 if(dabs(xn).gt.1.d-10) then 196 i1(1)=2 197 i1(2)=3 198 i2(1)=1 199 i2(2)=1 200 x1(1)=xn 201 x1(2)=xn 202 x2(1)=-yn 203 x2(2)=-zn 204 elseif(dabs(yn).gt.1.d-10) then 205 i1(1)=1 206 i1(2)=3 207 i2(1)=2 208 i2(2)=2 209 x1(1)=yn 210 x1(2)=yn 211 x2(1)=-xn 212 x2(2)=-zn 213 elseif(dabs(zn).gt.1.d-10) then 214 i1(1)=1 215 i1(2)=2 216 i2(1)=3 217 i2(2)=3 218 x1(1)=zn 219 x1(2)=zn 220 x2(1)=-xn 221 x2(2)=-yn 222 endif 223! 224! generating two MPC's expressing that the nodes cannot 225! move in planes perpendicular to the cyclic symmetry 226! axis 227! 228 do k=1,2 229 idof=8*(node-1)+i1(k) 230 call nident(ikmpc,idof,nmpc,id) 231 if(id.gt.0) then 232 if(ikmpc(id).eq.idof) then 233 write(*,*) 234 & '*ERROR reading *SELECT CYCLIC SYMMETRY MODES:' 235 write(*,*) ' node',node, 236 & ' on cyclic symmetry' 237 write(*,*) ' axis is used in other MPC' 238 ier=1 239 return 240 endif 241 endif 242 nmpc=nmpc+1 243 ipompc(nmpc)=mpcfree 244 labmpc(nmpc)=' ' 245! 246! updating ikmpc and ilmpc 247! 248 do j=nmpc,id+2,-1 249 ikmpc(j)=ikmpc(j-1) 250 ilmpc(j)=ilmpc(j-1) 251 enddo 252 ikmpc(id+1)=idof 253 ilmpc(id+1)=nmpc 254! 255 nodempc(1,mpcfree)=node 256 nodempc(2,mpcfree)=i1(k) 257 coefmpc(mpcfree)=x1(k) 258 mpcfree=nodempc(3,mpcfree) 259 if(mpcfree.eq.0) then 260 write(*,*) 261 & '*ERROR reading *SELECT CYCLIC SYMMETRY MODES:' 262 write(*,*) ' increase memmpc_' 263 ier=1 264 return 265 endif 266 nodempc(1,mpcfree)=node 267 nodempc(2,mpcfree)=i2(k) 268 coefmpc(mpcfree)=x2(k) 269 mpcfreeold=mpcfree 270 mpcfree=nodempc(3,mpcfree) 271 if(mpcfree.eq.0) then 272 write(*,*) 273 & '*ERROR reading *SELECT CYCLIC SYMMETRY MODES:' 274 write(*,*) ' increase memmpc_' 275 ier=1 276 return 277 endif 278 nodempc(3,mpcfreeold)=0 279 enddo 280 elseif(ns(2).eq.1) then 281! 282! nodal diameter 1 283! 284 if(dabs(xn).gt.1.d-10) then 285 i3=1 286 i4=2 287 i5=3 288 x3=xn 289 x4=yn 290 x5=zn 291 elseif(dabs(yn).gt.1.d-10) then 292 i3=2 293 i4=2 294 i5=3 295 x3=yn 296 x4=xn 297 x5=zn 298 else 299 i3=3 300 i4=1 301 i5=2 302 x3=zn 303 x4=xn 304 x5=yn 305 endif 306! 307! generating one MPC expressing that the nodes should 308! not move along the axis 309! 310 idof=8*(node-1)+i3 311 call nident(ikmpc,idof,nmpc,id) 312 if(id.gt.0) then 313 if(ikmpc(id).eq.idof) then 314 write(*,*) 315 & '*ERROR reading *SELECT CYCLIC SYMMETRY MODES:' 316 write(*,*) ' node',node, 317 & ' on cyclic symmetry' 318 write(*,*) ' axis is used in other MPC' 319 ier=1 320 return 321 endif 322 endif 323 nmpc=nmpc+1 324 ipompc(nmpc)=mpcfree 325 labmpc(nmpc)=' ' 326! 327! updating ikmpc and ilmpc 328! 329 do j=nmpc,id+2,-1 330 ikmpc(j)=ikmpc(j-1) 331 ilmpc(j)=ilmpc(j-1) 332 enddo 333 ikmpc(id+1)=idof 334 ilmpc(id+1)=nmpc 335! 336 nodempc(1,mpcfree)=node 337 nodempc(2,mpcfree)=i3 338 coefmpc(mpcfree)=x3 339 mpcfree=nodempc(3,mpcfree) 340 if(mpcfree.eq.0) then 341 write(*,*) 342 & '*ERROR reading *SELECT CYCLIC SYMMETRY MODES:' 343 write(*,*) ' increase memmpc_' 344 ier=1 345 return 346 endif 347 nodempc(1,mpcfree)=node 348 nodempc(2,mpcfree)=i4 349 coefmpc(mpcfree)=x4 350 mpcfree=nodempc(3,mpcfree) 351 if(mpcfree.eq.0) then 352 write(*,*) 353 & '*ERROR reading *SELECT CYCLIC SYMMETRY MODES:' 354 write(*,*) ' increase memmpc_' 355 ier=1 356 return 357 endif 358 nodempc(1,mpcfree)=node 359 nodempc(2,mpcfree)=i5 360 coefmpc(mpcfree)=x5 361 mpcfreeold=mpcfree 362 mpcfree=nodempc(3,mpcfree) 363 if(mpcfree.eq.0) then 364 write(*,*) 365 & '*ERROR reading *SELECT CYCLIC SYMMETRY MODES:' 366 write(*,*) ' increase memmpc_' 367 ier=1 368 return 369 endif 370 nodempc(3,mpcfreeold)=0 371 else 372 do k=1,3 373 idof=8*(node-1)+k 374 call nident(ikmpc,idof,nmpc,id) 375 if(id.gt.0) then 376 if(ikmpc(id).eq.idof) then 377 write(*,*) 378 & '*ERROR reading *SELECT CYCLIC SYMMETRY MODES:' 379 write(*,*) ' node',node, 380 & ' on cyclic symmetry' 381 write(*,*) ' axis is used in other MPC' 382 ier=1 383 return 384 endif 385 endif 386 nmpc=nmpc+1 387 ipompc(nmpc)=mpcfree 388 labmpc(nmpc)=' ' 389! 390! updating ikmpc and ilmpc 391! 392 do j=nmpc,id+2,-1 393 ikmpc(j)=ikmpc(j-1) 394 ilmpc(j)=ilmpc(j-1) 395 enddo 396 ikmpc(id+1)=idof 397 ilmpc(id+1)=nmpc 398! 399 nodempc(1,mpcfree)=node 400 nodempc(2,mpcfree)=k 401 coefmpc(mpcfree)=1.d0 402 mpcfreeold=mpcfree 403 mpcfree=nodempc(3,mpcfree) 404 if(mpcfree.eq.0) then 405 write(*,*) 406 & '*ERROR reading *SELECT CYCLIC SYMMETRY MODES:' 407 write(*,*) ' increase memmpc_' 408 ier=1 409 return 410 endif 411 nodempc(3,mpcfreeold)=0 412 enddo 413 endif 414 endif 415 enddo 416! 417 cs(2,ij)=ns(2)+0.5d0 418 cs(3,ij)=ns(3)+0.5d0 419 enddo 420! 421 if(irepeat.eq.0) irepeat=1 422! 423 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 424 & ipoinp,inp,ipoinpc) 425! 426c do j=1,nmpc 427c call writempc(ipompc,nodempc,coefmpc,labmpc,j) 428c enddo 429! 430 return 431 end 432 433