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 rhsnodef(co,kon,ne, 20 & ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc, 21 & nforc,fext,nactdof,nmethod,ikmpc,ntmat_,iperturb, 22 & mi,ikactmech,nactmech,ntrans,inotr,trab) 23! 24! filling the right hand side load vector b 25! 26! b contains the contributions due to mechanical forces only 27! 28 implicit none 29! 30 integer kon(*),ipompc(*),nodempc(3,*), 31 & nodeforc(2,*),ndirforc(*),ikmpc(*),mi(*), 32 & nactdof(0:mi(2),*),konl(20), 33 & idummy,rhsi,ntrans,inotr(2,*), 34 & ne,nmpc,nforc,nmethod,m,idm,idir,itr, 35 & iperturb(*),i,j,k,idist,jj,node, 36 & id,ist,index,jdof1,jdof,node1,ntmat_, 37 & icalccg,ikactmech(*),nactmech 38! 39 real*8 co(3,*),coefmpc(*),xforc(*), 40 & fext(*),a(3,3),fnext,trab(7,*) 41! 42 icalccg=0 43! 44! 45! point forces 46! 47! modal dynamics and steady state dynamics: 48! location of nonzeros is stored 49! 50 if((nmethod.ge.4).and.(iperturb(1).lt.2)) then 51 do i=1,nforc 52 if(ndirforc(i).gt.3) cycle 53 if(dabs(xforc(i)).lt.1.d-30) cycle 54! 55! checking for local coordinate system 56! 57 node=nodeforc(1,i) 58 if(ntrans.eq.0) then 59 itr=0 60 else 61 itr=inotr(1,node) 62 endif 63! 64 if(itr.eq.0) then 65! 66! no local coordinate system 67! 68 jdof=nactdof(ndirforc(i),node) 69 if(jdof.gt.0) then 70 fext(jdof)=fext(jdof)+xforc(i) 71 call nident(ikactmech,jdof-1,nactmech, 72 & idm) 73 do 74 if(idm.gt.0) then 75 if(ikactmech(idm).eq.jdof-1) exit 76 endif 77 nactmech=nactmech+1 78 do m=nactmech,idm+2,-1 79 ikactmech(m)=ikactmech(m-1) 80 enddo 81 ikactmech(idm+1)=jdof-1 82 exit 83 enddo 84 else 85! 86! node is a dependent node of a MPC: distribute 87! the forces among the independent nodes 88! (proportional to their coefficients) 89! 90 jdof=8*(node-1)+ndirforc(i) 91 call nident(ikmpc,jdof,nmpc,id) 92 if(id.gt.0) then 93 if(ikmpc(id).eq.jdof) then 94 ist=ipompc(id) 95 index=nodempc(3,ist) 96 if(index.eq.0) cycle 97 do 98 jdof=nactdof(nodempc(2,index), 99 & nodempc(1,index)) 100 if(jdof.gt.0) then 101 fext(jdof)=fext(jdof)- 102 & coefmpc(index)*xforc(i)/coefmpc(ist) 103 call nident(ikactmech,jdof-1,nactmech, 104 & idm) 105 do 106 if(idm.gt.0) then 107 if(ikactmech(idm).eq.jdof-1) exit 108 endif 109 nactmech=nactmech+1 110 do m=nactmech,idm+2,-1 111 ikactmech(m)=ikactmech(m-1) 112 enddo 113 ikactmech(idm+1)=jdof-1 114 exit 115 enddo 116 endif 117 index=nodempc(3,index) 118 if(index.eq.0) exit 119 enddo 120 endif 121 endif 122 endif 123 else 124! 125! local coordinate system 126! 127 call transformatrix(trab(1,itr),co(1,node),a) 128 idir=ndirforc(i) 129! 130! loop over all dofs 131! 132 do j=1,3 133 jdof=nactdof(j,node) 134 if(jdof.gt.0) then 135 fext(jdof)=fext(jdof)+a(j,idir)*xforc(i) 136 call nident(ikactmech,jdof-1,nactmech, 137 & idm) 138 do 139 if(idm.gt.0) then 140 if(ikactmech(idm).eq.jdof-1) exit 141 endif 142 nactmech=nactmech+1 143 do m=nactmech,idm+2,-1 144 ikactmech(m)=ikactmech(m-1) 145 enddo 146 ikactmech(idm+1)=jdof-1 147 exit 148 enddo 149 else 150! 151! node is a dependent node of a MPC: distribute 152! the forces among the independent nodes 153! (proportional to their coefficients) 154! 155 jdof=8*(node-1)+j 156 call nident(ikmpc,jdof,nmpc,id) 157 if(id.gt.0) then 158 if(ikmpc(id).eq.jdof) then 159 ist=ipompc(id) 160 index=nodempc(3,ist) 161 if(index.eq.0) cycle 162 do 163 jdof=nactdof(nodempc(2,index), 164 & nodempc(1,index)) 165 if(jdof.gt.0) then 166 fext(jdof)=fext(jdof)- 167 & coefmpc(index)*a(j,idir)*xforc(i)/ 168 & coefmpc(ist) 169 call nident(ikactmech,jdof-1,nactmech, 170 & idm) 171 do 172 if(idm.gt.0) then 173 if(ikactmech(idm).eq.jdof-1) exit 174 endif 175 nactmech=nactmech+1 176 do m=nactmech,idm+2,-1 177 ikactmech(m)=ikactmech(m-1) 178 enddo 179 ikactmech(idm+1)=jdof-1 180 exit 181 enddo 182 endif 183 index=nodempc(3,index) 184 if(index.eq.0) exit 185 enddo 186 endif 187 endif 188 endif 189 enddo 190 191 192 endif 193 enddo 194 else 195! 196! other procedures 197! 198 idummy=0 199 rhsi=1 200 call mafillsmforc(nforc,ndirforc,nodeforc,xforc,nactdof, 201 & fext,ipompc,nodempc,coefmpc,mi,rhsi,fnext,idummy, 202 & ntrans,inotr,trab,co) 203! 204c do i=1,nforc 205c if(ndirforc(i).gt.3) cycle 206c jdof=nactdof(ndirforc(i),nodeforc(1,i)) 207c if(jdof.gt.0) then 208c fext(jdof)=fext(jdof)+xforc(i) 209c else 210c! 211c! node is a dependent node of a MPC: distribute 212c! the forces among the independent nodes 213c! (proportional to their coefficients) 214c! 215c jdof=8*(nodeforc(1,i)-1)+ndirforc(i) 216c call nident(ikmpc,jdof,nmpc,id) 217c if(id.gt.0) then 218c if(ikmpc(id).eq.jdof) then 219c ist=ipompc(id) 220c index=nodempc(3,ist) 221c if(index.eq.0) cycle 222c do 223c jdof=nactdof(nodempc(2,index),nodempc(1,index)) 224c if(jdof.gt.0) then 225c fext(jdof)=fext(jdof)- 226c & coefmpc(index)*xforc(i)/coefmpc(ist) 227c endif 228c index=nodempc(3,index) 229c if(index.eq.0) exit 230c enddo 231c endif 232c endif 233c endif 234c enddo 235 endif 236! 237 return 238 end 239