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 shape8humass(xi,et,ze,xl,xsj,shp,iflag)
20!
21!     shape functions and derivatives for a 8-node linear isoparametric
22!     solid element
23!
24!     iflag=1: calculate only the value of the shape functions
25!     iflag=2: calculate the value of the shape functions and
26!              the Jacobian determinant
27!     iflag=3: calculate the value of the shape functions, the
28!              value of their derivatives w.r.t. the global
29!              coordinates and the Jacobian determinant
30!
31!     author: Otto-Ernst Bernhardi
32!
33      implicit none
34!
35      integer i,j,k,iflag
36!
37      real*8 shp(4,23),xs(3,3),xsi(3,3),xl(3,23),sh(3),xsi0(3,3)
38!
39      real*8 xi,et,ze,xsj,omg,omh,omr,opg,oph,opr
40!
41!
42!
43      if(iflag.gt.2) then
44!
45!        local derivatives at center point: xi-derivative
46!
47         shp(1,1)=-1.0d0/8.d0
48         shp(1,2)=1.0d0/8.d0
49         shp(1,3)=1.0d0/8.d0
50         shp(1,4)=-1.0d0/8.d0
51         shp(1,5)=-1.0d0/8.d0
52         shp(1,6)=1.0d0/8.d0
53         shp(1,7)=1.0d0/8.d0
54         shp(1,8)=-1.0d0/8.d0
55!
56!        local derivatives at center point: eta-derivative
57!
58         shp(2,1)=-1.0d0/8.d0
59         shp(2,2)=-1.0d0/8.d0
60         shp(2,3)=1.0d0/8.d0
61         shp(2,4)=1.0d0/8.d0
62         shp(2,5)=-1.0d0/8.d0
63         shp(2,6)=-1.0d0/8.d0
64         shp(2,7)=1.0d0/8.d0
65         shp(2,8)=1.0d0/8.d0
66!
67!        local derivatives at center point: zeta-derivative
68!
69         shp(3,1)=-1.0d0/8.d0
70         shp(3,2)=-1.0d0/8.d0
71         shp(3,3)=-1.0d0/8.d0
72         shp(3,4)=-1.0d0/8.d0
73         shp(3,5)=1.0d0/8.d0
74         shp(3,6)=1.0d0/8.d0
75         shp(3,7)=1.0d0/8.d0
76         shp(3,8)=1.0d0/8.d0
77!
78!        computation of the local derivative of the global coordinates
79!        (xs) at center point of element.
80!
81         do i=1,3
82            do j=1,3
83               xs(i,j)=0.d0
84               do k=1,8
85                  xs(i,j)=xs(i,j)+xl(i,k)*shp(j,k)
86               enddo
87            enddo
88         enddo
89!
90!        computation of the jacobian determinant at center point
91!
92         xsj=xs(1,1)*(xs(2,2)*xs(3,3)-xs(2,3)*xs(3,2))
93     &        -xs(1,2)*(xs(2,1)*xs(3,3)-xs(2,3)*xs(3,1))
94     &        +xs(1,3)*(xs(2,1)*xs(3,2)-xs(2,2)*xs(3,1))
95!
96!        computation of the global derivative of the local coordinates
97!        at center point of element.
98!
99         xsi0(1,1)=(xs(2,2)*xs(3,3)-xs(3,2)*xs(2,3))/xsj
100         xsi0(1,2)=(xs(1,3)*xs(3,2)-xs(1,2)*xs(3,3))/xsj
101         xsi0(1,3)=(xs(1,2)*xs(2,3)-xs(2,2)*xs(1,3))/xsj
102         xsi0(2,1)=(xs(2,3)*xs(3,1)-xs(2,1)*xs(3,3))/xsj
103         xsi0(2,2)=(xs(1,1)*xs(3,3)-xs(3,1)*xs(1,3))/xsj
104         xsi0(2,3)=(xs(1,3)*xs(2,1)-xs(1,1)*xs(2,3))/xsj
105         xsi0(3,1)=(xs(2,1)*xs(3,2)-xs(3,1)*xs(2,2))/xsj
106         xsi0(3,2)=(xs(1,2)*xs(3,1)-xs(1,1)*xs(3,2))/xsj
107         xsi0(3,3)=(xs(1,1)*xs(2,2)-xs(2,1)*xs(1,2))/xsj
108!
109      endif
110!
111!     shape functions and their global derivatives
112!
113      omg=1.d0-xi
114      omh=1.d0-et
115      omr=1.d0-ze
116      opg=1.d0+xi
117      oph=1.d0+et
118      opr=1.d0+ze
119!
120!     shape functions
121!
122      shp(4, 1)=omg*omh*omr/8.d0
123      shp(4, 2)=opg*omh*omr/8.d0
124      shp(4, 3)=opg*oph*omr/8.d0
125      shp(4, 4)=omg*oph*omr/8.d0
126      shp(4, 5)=omg*omh*opr/8.d0
127      shp(4, 6)=opg*omh*opr/8.d0
128      shp(4, 7)=opg*oph*opr/8.d0
129      shp(4, 8)=omg*oph*opr/8.d0
130c
131c     change on 190315: set shape functions to
132c     zero in order to obtain convergence in
133c     contact calculations with c3d8i
134c
135c      shp(4, 9)=0.d0
136c      shp(4,10)=0.d0
137c      shp(4,11)=0.d0
138      shp(4, 9)=1.0d0-xi*xi
139      shp(4,10)=1.0d0-et*et
140      shp(4,11)=1.0d0-ze*ze
141!
142      if(iflag.eq.1) return
143!
144!     local derivatives of the shape functions: xi-derivative
145!
146      shp(1, 1)=-omh*omr/8.d0
147      shp(1, 2)=omh*omr/8.d0
148      shp(1, 3)=oph*omr/8.d0
149      shp(1, 4)=-oph*omr/8.d0
150      shp(1, 5)=-omh*opr/8.d0
151      shp(1, 6)=omh*opr/8.d0
152      shp(1, 7)=oph*opr/8.d0
153      shp(1, 8)=-oph*opr/8.d0
154      shp(1, 9)=-2.d0*xi
155      shp(1,10)=0.d0
156      shp(1,11)=0.d0
157!
158!     local derivatives of the shape functions: eta-derivative
159!
160      shp(2, 1)=-omg*omr/8.d0
161      shp(2, 2)=-opg*omr/8.d0
162      shp(2, 3)=opg*omr/8.d0
163      shp(2, 4)=omg*omr/8.d0
164      shp(2, 5)=-omg*opr/8.d0
165      shp(2, 6)=-opg*opr/8.d0
166      shp(2, 7)=opg*opr/8.d0
167      shp(2, 8)=omg*opr/8.d0
168      shp(2, 9)=0.d0
169      shp(2,10)=-2.d0*et
170      shp(2,11)=0.d0
171!
172!     local derivatives of the shape functions: zeta-derivative
173!
174      shp(3, 1)=-omg*omh/8.d0
175      shp(3, 2)=-opg*omh/8.d0
176      shp(3, 3)=-opg*oph/8.d0
177      shp(3, 4)=-omg*oph/8.d0
178      shp(3, 5)=omg*omh/8.d0
179      shp(3, 6)=opg*omh/8.d0
180      shp(3, 7)=opg*oph/8.d0
181      shp(3, 8)=omg*oph/8.d0
182      shp(3, 9)=0.d0
183      shp(3,10)=0.d0
184      shp(3,11)=-2.d0*ze
185!
186!     computation of the local derivative of the global coordinates
187!     (xs). Incompatible modes are not included here.
188!
189      do i=1,3
190        do j=1,3
191          xs(i,j)=0.d0
192          do k=1,8
193            xs(i,j)=xs(i,j)+xl(i,k)*shp(j,k)
194          enddo
195        enddo
196      enddo
197!
198!     computation of the jacobian determinant
199!
200      xsj=xs(1,1)*(xs(2,2)*xs(3,3)-xs(2,3)*xs(3,2))
201     &   -xs(1,2)*(xs(2,1)*xs(3,3)-xs(2,3)*xs(3,1))
202     &   +xs(1,3)*(xs(2,1)*xs(3,2)-xs(2,2)*xs(3,1))
203!
204      if(iflag.eq.2) return
205!
206!     computation of the global derivative of the local coordinates
207!     (xsi) (inversion of xs)
208!
209      xsi(1,1)=(xs(2,2)*xs(3,3)-xs(3,2)*xs(2,3))/xsj
210      xsi(1,2)=(xs(1,3)*xs(3,2)-xs(1,2)*xs(3,3))/xsj
211      xsi(1,3)=(xs(1,2)*xs(2,3)-xs(2,2)*xs(1,3))/xsj
212      xsi(2,1)=(xs(2,3)*xs(3,1)-xs(2,1)*xs(3,3))/xsj
213      xsi(2,2)=(xs(1,1)*xs(3,3)-xs(3,1)*xs(1,3))/xsj
214      xsi(2,3)=(xs(1,3)*xs(2,1)-xs(1,1)*xs(2,3))/xsj
215      xsi(3,1)=(xs(2,1)*xs(3,2)-xs(3,1)*xs(2,2))/xsj
216      xsi(3,2)=(xs(1,2)*xs(3,1)-xs(1,1)*xs(3,2))/xsj
217      xsi(3,3)=(xs(1,1)*xs(2,2)-xs(2,1)*xs(1,2))/xsj
218!
219!     computation of the global derivatives of the shape functions
220!
221      do k=1,8
222        do j=1,3
223          sh(j)=shp(1,k)*xsi(1,j)+shp(2,k)*xsi(2,j)+shp(3,k)*xsi(3,j)
224        enddo
225        do j=1,3
226          shp(j,k)=sh(j)
227        enddo
228      enddo
229      do k=9,11
230        do j=1,3
231          sh(j)=shp(1,k)*xsi0(1,j)+shp(2,k)*xsi0(2,j)+shp(3,k)*xsi0(3,j)
232        enddo
233        do j=1,3
234          shp(j,k)=sh(j)
235        enddo
236      enddo
237!
238      return
239      end
240