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 hgstiffness(s,elas,a,gs)
20!
21!     hourglass control stiffness for 8-node solid mean strain element
22!
23!     Reference: Flanagan, D.P., Belytschko, T.; "Uniform  strain hexahedron
24!     and quadrilateral with orthogonal Hourglass control". Int. J. Num.
25!     Meth. Engg., Vol. 17, 679-706, 1981.
26!
27!     author: Otto-Ernst Bernhardi
28!
29      implicit none
30!
31      integer ii1,jj1,ii,jj,m1
32!
33      real*8 s(60,60),gs(8,4),a,elas(1),hgls,ahr
34!
35!
36!
37      ahr=elas(1)*a
38c     write(6,*) "stiffness:", ahr
39!
40      jj1=1
41      do jj=1,8
42         ii1=1
43         do ii=1,jj
44           hgls=0.0d0
45           do m1=1,4
46              hgls=hgls+gs(jj,m1)*gs(ii,m1)
47           enddo
48           hgls=hgls*ahr
49           s(ii1,jj1)=s(ii1,jj1)+hgls
50           s(ii1+1,jj1+1)=s(ii1+1,jj1+1)+hgls
51           s(ii1+2,jj1+2)=s(ii1+2,jj1+2)+hgls
52           ii1=ii1+3
53         enddo
54         jj1=jj1+3
55      enddo
56      return
57      end
58