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