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 calcspringforc(imat,elcon,nelcon,ncmat_,ntmat_,t1l, 20 & kode,plicon,nplicon,npmat_,senergy,nener,fk,val) 21! 22! calculates the spring forc and the spring energy (node-to-face penalty) 23! 24 implicit none 25! 26 integer i,imat,ncmat_,ntmat_,kode,niso,id,nplicon(0:ntmat_,*), 27 & npmat_,nelcon(2,*),nener 28! 29 real*8 t1l,elcon(0:ncmat_,ntmat_,*),elconloc(21),plconloc(802), 30 & xk,fk,val,xiso(200),yiso(200),plicon(0:2*npmat_,ntmat_,*), 31 & senergy 32! 33! 34! 35! interpolating the material data 36! 37 call materialdata_sp(elcon,nelcon,imat,ntmat_,i,t1l, 38 & elconloc,kode,plicon,nplicon,npmat_,plconloc,ncmat_) 39! 40! calculating the spring force and the spring constant 41! 42 if(kode.eq.2)then 43 xk=elconloc(1) 44 fk=xk*val 45 if(nener.eq.1) then 46 senergy=fk*val/2.d0 47 endif 48 else 49 niso=int(plconloc(801)) 50 do i=1,niso 51 xiso(i)=plconloc(2*i-1) 52 yiso(i)=plconloc(2*i) 53 enddo 54 call ident(xiso,val,niso,id) 55 if(id.eq.0) then 56 xk=0.d0 57 fk=yiso(1) 58 if(nener.eq.1) then 59 senergy=fk*val 60 endif 61 elseif(id.eq.niso) then 62 xk=0.d0 63 fk=yiso(niso) 64 if(nener.eq.1) then 65 senergy=yiso(1)*xiso(1) 66 do i=2,niso 67 senergy=senergy+(xiso(i)-xiso(i-1))* 68 & (yiso(i)+yiso(i-1))/2.d0 69 enddo 70 senergy=senergy+(val-xiso(niso))*yiso(niso) 71 endif 72 else 73 xk=(yiso(id+1)-yiso(id))/(xiso(id+1)-xiso(id)) 74 fk=yiso(id)+xk*(val-xiso(id)) 75 if(nener.eq.1) then 76 senergy=yiso(1)*xiso(1) 77 do i=2, id 78 senergy=senergy+(xiso(i)-xiso(i-1))* 79 & (yiso(i)+yiso(i-1))/2.d0 80 enddo 81 senergy=senergy+(val-xiso(id))*(fk+yiso(id))/2.d0 82 endif 83 endif 84 endif 85! 86 return 87 end 88 89