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 forcesolve(zc,nev,a,b,x,eiga,eigb,eigxx,iter,d, 20 & neq,z,istartnmd,iendnmd,nmd,cyclicsymmetry,neqact, 21 & igeneralizedforce) 22! 23! solves for the complex eigenfrequencies due to Coriolis 24! forces 25! 26 implicit none 27! 28 logical wantx 29! 30 integer nev,neq,iter(*),i,j,k,l,istartnmd(*),iendnmd(*),nmd, 31 & cyclicsymmetry,neqact,igeneralizedforce 32! 33 real*8 z(neq,*),d(*) 34! 35 complex*16 a(nev,*),b(nev,*),x(nev,*),eiga(*),eigb(*),eigxx(*), 36 & zc(neqact,*) 37! 38 if(igeneralizedforce.eq.0) then 39! 40! no generalized force: multiplication with the eigenmodes 41! is necessary 42! 43 if(cyclicsymmetry.eq.0) then 44 do i=1,nev 45 do j=1,nev 46 do k=1,neq 47 a(i,j)=a(i,j)+z(k,i)*zc(k,j) 48 enddo 49 enddo 50 write(*,*) 51 & 'aerodynamic stiffness/structural stiffness = ', 52 & a(i,i)/d(i) 53 a(i,i)=a(i,i)+d(i) 54 b(i,i)=(1.d0,0.d0) 55 enddo 56 else 57! 58! cyclic symmetry 59! 60 do l=1,nmd 61 do i=istartnmd(l),iendnmd(l) 62 do j=istartnmd(l),iendnmd(l) 63 do k=1,neqact 64 a(i,j)=a(i,j)+z(k,i)*zc(k,j)- 65 & z(k+neqact,i)*zc(k,j)*(0.d0,1.d0) 66 enddo 67 enddo 68 write(*,*) 69 & 'aerodynamic stiffness/structural stiffness = ', 70 & a(i,i)/d(i) 71 a(i,i)=a(i,i)+d(i) 72 b(i,i)=(1.d0,0.d0) 73 enddo 74 enddo 75 endif 76 else 77! 78! generalized force: the a-matrix is (apart from the diagonal) 79! known 80! 81 if(cyclicsymmetry.eq.0) then 82 do i=1,nev 83 write(*,*) 84 & 'aerodynamic stiffness/structural stiffness = ', 85 & a(i,i)/d(i) 86 a(i,i)=a(i,i)+d(i) 87 b(i,i)=(1.d0,0.d0) 88 enddo 89 else 90! 91! cyclic symmetry 92! 93 do l=1,nmd 94 do i=istartnmd(l),iendnmd(l) 95 write(*,*) 96 & 'aerodynamic stiffness/structural stiffness = ', 97 & a(i,i)/d(i) 98 a(i,i)=a(i,i)+d(i) 99 b(i,i)=(1.d0,0.d0) 100 enddo 101 enddo 102 endif 103 endif 104! 105 wantx=.true. 106! 107! solving for the complex eigenvalues 108! 109 call dlzhes(nev,a,nev,b,nev,x,nev,wantx) 110 call dlzit(nev,a,nev,b,nev,x,nev,wantx,iter,eiga,eigb) 111! 112 do i=1,nev 113 if(iter(i).eq.-1) then 114 write(*,*) '*ERROR in forcesolve: fatal error' 115 write(*,*) ' in dlzit' 116 call exit(201) 117 elseif(cdabs(eigb(i)).lt.1.d-10) then 118 write(*,*) '*ERROR in forcesolve: eigenvalue' 119 write(*,*) ' out of bounds' 120 call exit(201) 121 else 122 eigxx(i)=cdsqrt(eiga(i)/eigb(i)) 123 endif 124 enddo 125! 126 return 127 end 128 129 130