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