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 modaldynamics(inpc,textpart,nmethod,tinc,tper,iexpl, 20 & istep,istat,n,iline,ipol,inl,ipoinp,inp,iperturb,isolver, 21 & cs,mcs,ipoinpc,idrct,ctrl,tmin,tmax,nforc,nload,nbody,iprestr, 22 & t0,t1,ithermal,nk,vold,veold,xmodal,set,nset,mi,cyclicsymmetry, 23 & ier) 24! 25! reading the input deck: *MODAL DYNAMIC 26! 27 implicit none 28! 29 logical steadystate,nodalset 30! 31 character*1 inpc(*) 32 character*20 solver 33 character*81 set(*),noset 34 character*132 textpart(16) 35! 36 integer nmethod,istep,istat,n,key,iexpl,iline,ipol,inl, 37 & ipoinp(2,*),inp(3,*),iperturb(*),isolver,i,mcs,ipoinpc(0:*), 38 & idrct,nforc,nload,nbody,iprestr,ithermal(*),j,nk,ipos,nset, 39 & mi(*),cyclicsymmetry,ier 40! 41 real*8 tinc,tper,cs(17,*),ctrl(*),tmin,tmax,t0(*),t1(*), 42 & vold(0:mi(2),*),veold(0:mi(2),*),xmodal(*) 43! 44 iexpl=0 45c iperturb(1)=0 46 iperturb(2)=0 47 idrct=1 48 tmin=0.d0 49 tmax=0.d0 50 steadystate=.false. 51 if((mcs.ne.0).and.(cs(2,1).ge.0.d0)) then 52 cyclicsymmetry=1 53 endif 54 nodalset=.false. 55! 56 if(istep.lt.1) then 57 write(*,*) 58 & '*ERROR reading *MODAL DYNAMIC: *MODAL DYNAMIC can only' 59 write(*,*) ' be used within a STEP' 60 ier=1 61 return 62 endif 63! 64! default solver 65! 66 solver=' ' 67 if(isolver.eq.0) then 68 solver(1:7)='SPOOLES' 69 elseif(isolver.eq.2) then 70 solver(1:16)='ITERATIVESCALING' 71 elseif(isolver.eq.3) then 72 solver(1:17)='ITERATIVECHOLESKY' 73 elseif(isolver.eq.4) then 74 solver(1:3)='SGI' 75 elseif(isolver.eq.5) then 76 solver(1:5)='TAUCS' 77 elseif(isolver.eq.7) then 78 solver(1:7)='PARDISO' 79 elseif(isolver.eq.8) then 80 solver(1:6)='PASTIX' 81 endif 82! 83 do i=2,n 84 if(textpart(i)(1:7).eq.'SOLVER=') then 85 read(textpart(i)(8:27),'(a20)') solver 86 elseif(textpart(i)(1:9).eq.'DIRECT=NO') then 87 idrct=0 88 elseif(textpart(i)(1:7).eq.'DELTMX=') then 89 read(textpart(i)(8:27),'(f20.0)',iostat=istat) ctrl(27) 90 elseif(textpart(i)(1:11).eq.'STEADYSTATE') then 91 steadystate=.true. 92 else 93 write(*,*) 94 & '*WARNING reading *MODAL DYNAMIC: parameter not recognized:' 95 write(*,*) ' ', 96 & textpart(i)(1:index(textpart(i),' ')-1) 97 call inputwarning(inpc,ipoinpc,iline, 98 &"*MODAL DYNAMIC%") 99 endif 100 enddo 101! 102 if(solver(1:7).eq.'SPOOLES') then 103 isolver=0 104 elseif(solver(1:16).eq.'ITERATIVESCALING') then 105 write(*,*) 106 & '*WARNING reading *MODAL DYNAMIC: the iterative scaling' 107 write(*,*) ' procedure is not available for modal' 108 write(*,*) ' dynamic calculations; the default solver' 109 write(*,*) ' is used' 110 elseif(solver(1:17).eq.'ITERATIVECHOLESKY') then 111 write(*,*) 112 & '*WARNING reading *MODAL DYNAMIC: the iterative scaling' 113 write(*,*) ' procedure is not available for modal' 114 write(*,*) ' dynamic calculations; the default solver' 115 write(*,*) ' is used' 116 elseif(solver(1:3).eq.'SGI') then 117 isolver=4 118 elseif(solver(1:5).eq.'TAUCS') then 119 isolver=5 120 elseif(solver(1:7).eq.'PARDISO') then 121 isolver=7 122 elseif(solver(1:6).eq.'PARDISO') then 123 isolver=8 124 else 125 write(*,*) '*WARNING reading *MODAL DYNAMIC: unknown solver;' 126 write(*,*) ' the default solver is used' 127 endif 128! 129c if((isolver.eq.2).or.(isolver.eq.3)) then 130c write(*,*) '*ERROR reading *MODAL DYNAMIC: the default solver ', 131c & solver 132c write(*,*) ' cannot be used for modal dynamic' 133c write(*,*) ' calculations ' 134c ier=1 135c return 136c endif 137! 138 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 139 & ipoinp,inp,ipoinpc) 140 if((istat.lt.0).or.(key.eq.1)) then 141 write(*,*) 142 & '*ERROR reading *MODAL DYNAMIC: definition not complete' 143 write(*,*) ' ' 144 call inputerror(inpc,ipoinpc,iline, 145 & "*MODAL DYNAMIC%",ier) 146 return 147 endif 148 read(textpart(1)(1:20),'(f20.0)',iostat=istat)tinc 149 if(istat.gt.0) then 150 call inputerror(inpc,ipoinpc,iline, 151 & "*MODAL DYNAMIC%",ier) 152 return 153 endif 154 read(textpart(2)(1:20),'(f20.0)',iostat=istat)tper 155 if(istat.gt.0) then 156 call inputerror(inpc,ipoinpc,iline, 157 & "*MODAL DYNAMIC%",ier) 158 return 159 endif 160 read(textpart(3)(1:20),'(f20.0)',iostat=istat) tmin 161 if(istat.gt.0) then 162 call inputerror(inpc,ipoinpc,iline, 163 & "*MODAL DYNAMIC%",ier) 164 return 165 endif 166 read(textpart(4)(1:20),'(f20.0)',iostat=istat) tmax 167 if(istat.gt.0) then 168 call inputerror(inpc,ipoinpc,iline, 169 & "*MODAL DYNAMIC%",ier) 170 return 171 endif 172! 173 if(steadystate) then 174! 175! modal dynamics calculation till steady state 176! 177 if(tper.le.0.d0) then 178 write(*,*) '*ERROR reading *MODAL DYNAMIC: relative error' 179 write(*,*) ' is nonpositive' 180 ier=1 181 return 182 endif 183 tper=-tper 184 if(tinc.le.0.d0) then 185 write(*,*) 186 & '*ERROR reading *MODAL DYNAMIC: initial increment' 187 write(*,*) ' size is nonpositive' 188 ier=1 189 return 190 endif 191 if(tmin.lt.0.d0) then 192 tmin=1.d-10 193 endif 194 if(tmax.lt.1.d-10) then 195 tmax=1.d+30 196 endif 197 else 198! 199! transient modal dynamics calculation 200! 201 if(tper.lt.0.d0) then 202 write(*,*) 203 & '*ERROR reading *MODAL DYNAMIC: step size is negative' 204 ier=1 205 return 206 elseif(tper.le.0.d0) then 207 tper=1.d0 208 endif 209 if(tinc.lt.0.d0) then 210 write(*,*) 211 & '*ERROR reading *MODAL DYNAMIC: initial increment size 212 &is negative' 213 ier=1 214 return 215 elseif(tinc.le.0.d0) then 216 tinc=tper 217 endif 218 if(tinc.gt.tper) then 219 write(*,*) 220 & '*ERROR reading *MODAL DYNAMIC: initial increment size 221 &exceeds step size' 222 ier=1 223 return 224 endif 225! 226 if(idrct.ne.1) then 227 if(tmin.lt.1.d-10*tper) then 228 tmin=min(tinc,1.d-10*tper) 229 endif 230 if(tmax.lt.1.d-10) then 231 tmax=1.d+30 232 endif 233 endif 234 endif 235! 236! removing the present loading 237! 238c nforc=0 239c nload=0 240c nbody=0 241c iprestr=0 242c if((ithermal.eq.1).or.(ithermal.eq.3)) then 243c do j=1,nk 244c t1(j)=t0(j) 245c enddo 246c endif 247! 248! resetting fields vold and veold after a frequency or 249! buckling step 250! 251 if((nmethod.eq.2).or.(nmethod.eq.3)) then 252 do i=1,nk 253 do j=1,3 254 vold(j,i)=0.d0 255 veold(j,i)=0.d0 256 enddo 257 enddo 258 endif 259! 260 nmethod=4 261! 262! correction for cyclic symmetric structures: 263! if the present step was not preceded by a frequency step 264! no nodal diameter has been selected. To make sure that 265! mastructcs is called instead of mastruct a fictitious 266! minimum nodal diameter is stored 267! 268 if((cyclicsymmetry.eq.1).and.(mcs.ne.0).and.(cs(2,1).lt.0.d0)) 269 & cs(2,1)=0.d0 270! 271 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 272 & ipoinp,inp,ipoinpc) 273! 274 return 275 end 276