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 subspace(d,aa,bb,cc,alpham,betam,nev,xini, 20 & cd,cv,time,rwork,lrw,m,jout,rpar,bj,iwork,liw,iddebdf,bjp) 21! 22! solves the linear dynamic equations mapped on the subspace 23! of the eigenvectors (only if there are dashpots in the 24! model) 25! 26 implicit none 27! 28 save time0 29! 30 integer nev,nev2,info(15),idid,lrw,iwork(*),liw,jout,id, 31 & iab,iaa,ibb,i,j,m,iddebdf 32! 33 real*8 d(*),aa(nev),bb(nev),cc(nev,*),alpham,betam, 34 & xini(*),cd(*),cv(*),time,time0,rtol,atol,rwork(*),rpar(*), 35 & bj(*),bjp(*) 36! 37 external df,djac 38! 39! 40 nev2=2*nev 41! 42! transferring fields into global field rpar 43! (needed for subroutine fd) 44! rpar contains (field, size): m+0.5, 1 45! alpham, 1 46! betam, 1 47! cc, nev**2 48! d, nev 49! time 50! aa(1)..aa(nev), nev 51! bb(1)..bb(nev), nev 52! 53 if(m.eq.1) then 54 rpar(2)=alpham 55 rpar(3)=betam 56 do j=1,nev 57 do i=1,nev 58 rpar(3+(j-1)*nev+i)=cc(i,j) 59 enddo 60 enddo 61 id=3+nev*nev 62 do i=1,nev 63 rpar(id+i)=d(i) 64 enddo 65! 66! copying the initial conditions for the system of first order 67! differential equations 68! 69 do i=1,nev 70 xini(i)=cd(i) 71 xini(nev+i)=cv(i) 72 enddo 73 endif 74! 75 iaa=3+nev*(1+nev)+1 76 rpar(iaa)=time 77 ibb=iaa+nev 78 do i=1,nev 79 rpar(iaa+i)=aa(i) 80 rpar(ibb+i)=bb(i) 81 enddo 82! 83 do i=1,3 84 info(i)=0 85 enddo 86 info(4)=1 87 info(5)=1 88 info(6)=0 89 rwork(1)=time 90! 91! absolute and relative tolerance for dderkf 92! 93 rtol=1.d-5 94 atol=1.d-3 95! 96 if(iddebdf.eq.0) then 97 call ddeabm(df,nev2,time0,xini,time,info,rtol,atol,idid,rwork, 98 & lrw,iwork,liw,rpar,nev) 99! 100 if((idid.ne.2).and.(idid.ne.3)) then 101 write(*,*) 102 & '*WARNING in subspace: ddeabm did not converge properly' 103 write(*,*) ' idid= ',idid 104 write(*,*) ' switch to routine ddebdf' 105 iddebdf=2 106 return 107 endif 108 else 109 call ddebdf(df,nev2,time0,xini,time,info,rtol,atol,idid,rwork, 110 & lrw,iwork,liw,rpar,nev,djac) 111 if((idid.ne.2).and.(idid.ne.3)) then 112 write(*,*) 113 & '*ERROR in subspace: ddebdf did not converge properly' 114 write(*,*) ' idid= ',idid 115 call exit(201) 116 endif 117 endif 118! 119! copying the solution into field bj 120! 121 do i=1,nev 122 bj(i)=xini(i) 123 bjp(i)=xini(nev+i) 124 enddo 125! 126 return 127 end 128! 129! subroutine df expressing the first order derivative as a function 130! of time and the function itself 131! 132 subroutine df(x,u,uprime,rpar,nev) 133! 134 implicit none 135! 136 integer nev,i,j,id,iab,iaa,ibb,nev2p1,m 137! 138 real*8 rpar(*),x,u(*),uprime(*) 139! 140 iaa=4+nev*(nev+1) 141 ibb=iaa+nev 142 id=3+nev*nev 143! 144 do i=1,nev 145 uprime(i)=u(nev+i) 146 uprime(nev+i)=rpar(iaa+i)+x*rpar(ibb+i) 147 & -rpar(id+i)*rpar(id+i)*u(i) 148 & -(rpar(2)+rpar(3)*rpar(id+i)*rpar(id+i))*u(nev+i) 149! 150! contribution of the dashpots 151! 152 do j=1,nev 153 uprime(nev+i)=uprime(nev+i)-rpar(3+(j-1)*nev+i)*u(nev+j) 154 enddo 155 enddo 156! 157 return 158 end 159! 160! subroutine djac 161! 162 subroutine djac(x,u,pd,nrowpd,rpar,nev) 163! 164 implicit none 165! 166 integer nrowpd,nev,id,i,j 167! 168 real*8 rpar(*),x,u(*),pd(nrowpd,*) 169! 170 id=3+nev*nev 171! 172 do i=1,nev 173 pd(i,nev+i)=1.d0 174 pd(nev+i,i)=-rpar(id+i)*rpar(id+i) 175 pd(nev+i,nev+i)=-(rpar(2)+rpar(3)*rpar(id+i)*rpar(id+i)) 176! 177! contribution of the dashpots 178! 179 do j=1,nev 180 pd(nev+i,nev+j)=pd(nev+i,nev+j)-rpar(3+(j-1)*nev+i) 181 enddo 182 enddo 183! 184 return 185 end 186