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 shape20h(xi,et,ze,xl,xsj,shp,iflag)
20!
21!     shape functions and derivatives for a 20-node quadratic
22!     isoparametric brick element. -1<=xi,et,ze<=1
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      implicit none
32!
33      integer i,j,k,iflag
34!
35      real*8 shp(4,20),xs(3,3),xsi(3,3),xl(3,20),shpe(4,20),dd,
36     &  dd1,dd2,dd3
37!
38      real*8 xi,et,ze,xsj,omg,omh,omr,opg,oph,opr,
39     &  tpgphpr,tmgphpr,tmgmhpr,tpgmhpr,tpgphmr,tmgphmr,tmgmhmr,tpgmhmr,
40     &  omgopg,omhoph,omropr,omgmopg,omhmoph,omrmopr
41!
42!
43!
44!     shape functions and their glocal derivatives
45!
46      omg=1.d0-xi
47      omh=1.d0-et
48      omr=1.d0-ze
49      opg=1.d0+xi
50      oph=1.d0+et
51      opr=1.d0+ze
52      tpgphpr=opg+oph+ze
53      tmgphpr=omg+oph+ze
54      tmgmhpr=omg+omh+ze
55      tpgmhpr=opg+omh+ze
56      tpgphmr=opg+oph-ze
57      tmgphmr=omg+oph-ze
58      tmgmhmr=omg+omh-ze
59      tpgmhmr=opg+omh-ze
60      omgopg=omg*opg/4.d0
61      omhoph=omh*oph/4.d0
62      omropr=omr*opr/4.d0
63      omgmopg=(omg-opg)/4.d0
64      omhmoph=(omh-oph)/4.d0
65      omrmopr=(omr-opr)/4.d0
66!
67!     shape functions
68!
69      shp(4, 1)=-omg*omh*omr*tpgphpr/8.d0
70      shp(4, 2)=-opg*omh*omr*tmgphpr/8.d0
71      shp(4, 3)=-opg*oph*omr*tmgmhpr/8.d0
72      shp(4, 4)=-omg*oph*omr*tpgmhpr/8.d0
73      shp(4, 5)=-omg*omh*opr*tpgphmr/8.d0
74      shp(4, 6)=-opg*omh*opr*tmgphmr/8.d0
75      shp(4, 7)=-opg*oph*opr*tmgmhmr/8.d0
76      shp(4, 8)=-omg*oph*opr*tpgmhmr/8.d0
77      shp(4, 9)=omgopg*omh*omr
78      shp(4,10)=omhoph*opg*omr
79      shp(4,11)=omgopg*oph*omr
80      shp(4,12)=omhoph*omg*omr
81      shp(4,13)=omgopg*omh*opr
82      shp(4,14)=omhoph*opg*opr
83      shp(4,15)=omgopg*oph*opr
84      shp(4,16)=omhoph*omg*opr
85      shp(4,17)=omropr*omg*omh
86      shp(4,18)=omropr*opg*omh
87      shp(4,19)=omropr*opg*oph
88      shp(4,20)=omropr*omg*oph
89!
90      if(iflag.eq.1) return
91!
92!     local derivatives of the shape functions: xi-derivative
93!
94      shpe(1, 1)=omh*omr*(tpgphpr-omg)/8.d0
95      shpe(1, 2)=(opg-tmgphpr)*omh*omr/8.d0
96      shpe(1, 3)=(opg-tmgmhpr)*oph*omr/8.d0
97      shpe(1, 4)=oph*omr*(tpgmhpr-omg)/8.d0
98      shpe(1, 5)=omh*opr*(tpgphmr-omg)/8.d0
99      shpe(1, 6)=(opg-tmgphmr)*omh*opr/8.d0
100      shpe(1, 7)=(opg-tmgmhmr)*oph*opr/8.d0
101      shpe(1, 8)=oph*opr*(tpgmhmr-omg)/8.d0
102      shpe(1, 9)=omgmopg*omh*omr
103      shpe(1,10)=omhoph*omr
104      shpe(1,11)=omgmopg*oph*omr
105      shpe(1,12)=-omhoph*omr
106      shpe(1,13)=omgmopg*omh*opr
107      shpe(1,14)=omhoph*opr
108      shpe(1,15)=omgmopg*oph*opr
109      shpe(1,16)=-omhoph*opr
110      shpe(1,17)=-omropr*omh
111      shpe(1,18)=omropr*omh
112      shpe(1,19)=omropr*oph
113      shpe(1,20)=-omropr*oph
114!
115!     local derivatives of the shape functions: eta-derivative
116!
117      shpe(2, 1)=omg*omr*(tpgphpr-omh)/8.d0
118      shpe(2, 2)=opg*omr*(tmgphpr-omh)/8.d0
119      shpe(2, 3)=opg*(oph-tmgmhpr)*omr/8.d0
120      shpe(2, 4)=omg*(oph-tpgmhpr)*omr/8.d0
121      shpe(2, 5)=omg*opr*(tpgphmr-omh)/8.d0
122      shpe(2, 6)=opg*opr*(tmgphmr-omh)/8.d0
123      shpe(2, 7)=opg*(oph-tmgmhmr)*opr/8.d0
124      shpe(2, 8)=omg*(oph-tpgmhmr)*opr/8.d0
125      shpe(2, 9)=-omgopg*omr
126      shpe(2,10)=omhmoph*opg*omr
127      shpe(2,11)=omgopg*omr
128      shpe(2,12)=omhmoph*omg*omr
129      shpe(2,13)=-omgopg*opr
130      shpe(2,14)=omhmoph*opg*opr
131      shpe(2,15)=omgopg*opr
132      shpe(2,16)=omhmoph*omg*opr
133      shpe(2,17)=-omropr*omg
134      shpe(2,18)=-omropr*opg
135      shpe(2,19)=omropr*opg
136      shpe(2,20)=omropr*omg
137!
138!     local derivatives of the shape functions: zeta-derivative
139!
140      shpe(3, 1)=omg*omh*(tpgphpr-omr)/8.d0
141      shpe(3, 2)=opg*omh*(tmgphpr-omr)/8.d0
142      shpe(3, 3)=opg*oph*(tmgmhpr-omr)/8.d0
143      shpe(3, 4)=omg*oph*(tpgmhpr-omr)/8.d0
144      shpe(3, 5)=omg*omh*(opr-tpgphmr)/8.d0
145      shpe(3, 6)=opg*omh*(opr-tmgphmr)/8.d0
146      shpe(3, 7)=opg*oph*(opr-tmgmhmr)/8.d0
147      shpe(3, 8)=omg*oph*(opr-tpgmhmr)/8.d0
148      shpe(3, 9)=-omgopg*omh
149      shpe(3,10)=-omhoph*opg
150      shpe(3,11)=-omgopg*oph
151      shpe(3,12)=-omhoph*omg
152      shpe(3,13)=omgopg*omh
153      shpe(3,14)=omhoph*opg
154      shpe(3,15)=omgopg*oph
155      shpe(3,16)=omhoph*omg
156      shpe(3,17)=omrmopr*omg*omh
157      shpe(3,18)=omrmopr*opg*omh
158      shpe(3,19)=omrmopr*opg*oph
159      shpe(3,20)=omrmopr*omg*oph
160!
161!     computation of the local derivative of the global coordinates
162!     (xs)
163!
164      do i=1,3
165        do j=1,3
166          xs(i,j)=0.d0
167          do k=1,20
168            xs(i,j)=xs(i,j)+xl(i,k)*shpe(j,k)
169          enddo
170        enddo
171      enddo
172!
173!     computation of the jacobian determinant
174!
175      dd1=xs(2,2)*xs(3,3)-xs(2,3)*xs(3,2)
176      dd2=xs(2,3)*xs(3,1)-xs(2,1)*xs(3,3)
177      dd3=xs(2,1)*xs(3,2)-xs(2,2)*xs(3,1)
178      xsj=xs(1,1)*dd1+xs(1,2)*dd2+xs(1,3)*dd3
179!
180      if(iflag.eq.2) return
181!
182      dd=1.d0/xsj
183!
184!     computation of the global derivative of the local coordinates
185!     (xsi) (inversion of xs)
186!
187      xsi(1,1)=dd1*dd
188      xsi(1,2)=(xs(1,3)*xs(3,2)-xs(1,2)*xs(3,3))*dd
189      xsi(1,3)=(xs(1,2)*xs(2,3)-xs(2,2)*xs(1,3))*dd
190      xsi(2,1)=dd2*dd
191      xsi(2,2)=(xs(1,1)*xs(3,3)-xs(3,1)*xs(1,3))*dd
192      xsi(2,3)=(xs(1,3)*xs(2,1)-xs(1,1)*xs(2,3))*dd
193      xsi(3,1)=dd3*dd
194      xsi(3,2)=(xs(1,2)*xs(3,1)-xs(1,1)*xs(3,2))*dd
195      xsi(3,3)=(xs(1,1)*xs(2,2)-xs(2,1)*xs(1,2))*dd
196!
197!     computation of the global derivatives of the shape functions
198!
199      do k=1,20
200        do j=1,3
201          shp(j,k)=shpe(1,k)*xsi(1,j)+shpe(2,k)*xsi(2,j)
202     &          +shpe(3,k)*xsi(3,j)
203        enddo
204      enddo
205c      do k=1,20
206c        do j=1,3
207c          shp(j,k)=shpe(j,k)
208c        enddo
209c      enddo
210!
211      return
212      end
213