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 shape15w(xi,et,ze,xl,xsj,shp,iflag)
20!
21!     shape functions and derivatives for a 15-node quadratic
22!     isoparametric wedge element. 0<=xi,et<=1,-1<=ze<=1,xi+et<=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!
32!     Copyright (c) 2003 WB
33!
34!     Written February 2003 on the basis of the Guido's shape function files
35!
36      implicit none
37!
38      integer i,j,k,iflag
39!
40      real*8 shp(4,15),xs(3,3),xsi(3,3),xl(3,15),sh(3)
41!
42      real*8 xi,et,ze,xsj,a
43!
44!
45!
46!     shape functions and their glocal derivatives
47!
48      a=1.d0-xi-et
49!
50!     shape functions
51!
52      shp(4,1)=-0.5d0*a*(1.d0-ze)*(2.d0*xi+2.d0*et+ze)
53      shp(4,2)=0.5d0*xi*(1.d0-ze)*(2.d0*xi-2.d0-ze)
54      shp(4,3)=0.5d0*et*(1.d0-ze)*(2.d0*et-2.d0-ze)
55      shp(4,4)=-0.5d0*a*(1.d0+ze)*(2.d0*xi+2.d0*et-ze)
56      shp(4,5)=0.5d0*xi*(1.d0+ze)*(2.d0*xi-2.d0+ze)
57      shp(4,6)=0.5d0*et*(1.d0+ze)*(2.d0*et-2.d0+ze)
58      shp(4,7)=2.d0*xi*a*(1.d0-ze)
59      shp(4,8)=2.d0*xi*et*(1.d0-ze)
60      shp(4,9)=2.d0*et*a*(1.d0-ze)
61      shp(4,10)=2.d0*xi*a*(1.d0+ze)
62      shp(4,11)=2.d0*xi*et*(1.d0+ze)
63      shp(4,12)=2.d0*et*a*(1.d0+ze)
64      shp(4,13)= a*(1.d0-ze*ze)
65      shp(4,14)=xi*(1.d0-ze*ze)
66      shp(4,15)=et*(1.d0-ze*ze)
67!
68      if(iflag.eq.1) return
69!
70!     local derivatives of the shape functions: xi-derivative
71!
72      shp(1,1)=  0.5d0*(1.d0-ze)*(4.0d0*xi+4.0d0*et+ze-2.d0)
73      shp(1,2)=  0.5d0*(1.d0-ze)*(4.0d0*xi-ze-2.d0)
74      shp(1,3)=  0.d0
75      shp(1,4)=  0.5d0*(1.d0+ze)*(4.0d0*xi+4.0d0*et-ze-2.d0)
76      shp(1,5)=  0.5d0*(1.d0+ze)*(4.0d0*xi+ze-2.d0)
77      shp(1,6)=  0.d0
78      shp(1,7)=  2.d0*(1.d0-ze)*(1.d0-2.d0*xi-et)
79      shp(1,8)=  2.d0*et*(1.d0-ze)
80      shp(1,9)= -2.d0*et*(1.d0-ze)
81      shp(1,10)= 2.d0*(1.d0+ze)*(1.d0-2.d0*xi-et)
82      shp(1,11)= 2.d0*et*(1.d0+ze)
83      shp(1,12)= -2.d0*et*(1.d0+ze)
84      shp(1,13)= -(1.d0-ze*ze)
85      shp(1,14)=  (1.d0-ze*ze)
86      shp(1,15)=  0.d0
87!
88!     local derivatives of the shape functions: eta-derivative
89!
90      shp(2,1)=  0.5d0*(1.d0-ze)*(4.0d0*xi+4.0d0*et+ze-2.d0)
91      shp(2,2)= 0.d0
92      shp(2,3)= 0.5d0*(1.d0-ze)*(4.0d0*et-ze-2.d0)
93      shp(2,4)= 0.5d0*(1.d0+ze)*(4.0d0*xi+4.0d0*et-ze-2.d0)
94      shp(2,5)= 0.d0
95      shp(2,6)= 0.5d0*(1.d0+ze)*(4.0d0*et+ze-2.d0)
96      shp(2,7)=-2.d0*xi*(1.d0-ze)
97      shp(2,8)= 2.d0*xi*(1.d0-ze)
98      shp(2,9)= 2.d0*(1.d0-ze)*(1.d0-xi-2.d0*et)
99      shp(2,10)=-2.d0*xi*(1.d0+ze)
100      shp(2,11)= 2.d0*xi*(1.d0+ze)
101      shp(2,12)= 2.d0*(1.d0+ze)*(1.d0-xi-2.d0*et)
102      shp(2,13)=-(1.d0-ze*ze)
103      shp(2,14)= 0.0d0
104      shp(2,15)= (1.d0-ze*ze)
105!
106!     local derivatives of the shape functions: zeta-derivative
107!
108      shp(3,1)=  a*(xi+et+ze-0.5d0)
109      shp(3,2)= xi*(-xi+ze+0.5d0)
110      shp(3,3)= et*(-et+ze+0.5d0)
111      shp(3,4)=  a*(-xi-et+ze+0.5d0)
112      shp(3,5)= xi*(xi+ze-0.5d0)
113      shp(3,6)= et*(et+ze-0.5d0)
114      shp(3,7)= -2*xi*a
115      shp(3,8)= -2*xi*et
116      shp(3,9)= -2*et*a
117      shp(3,10)= 2*xi*a
118      shp(3,11)= 2*xi*et
119      shp(3,12)= 2*et*a
120      shp(3,13)=-2*a*ze
121      shp(3,14)=-2*xi*ze
122      shp(3,15)=-2*et*ze
123!
124!     computation of the local derivative of the global coordinates
125!     (xs)
126!
127      do i=1,3
128        do j=1,3
129          xs(i,j)=0.d0
130          do k=1,15
131            xs(i,j)=xs(i,j)+xl(i,k)*shp(j,k)
132          enddo
133        enddo
134      enddo
135!
136!     computation of the jacobian determinant
137!
138      xsj=xs(1,1)*(xs(2,2)*xs(3,3)-xs(2,3)*xs(3,2))
139     &   -xs(1,2)*(xs(2,1)*xs(3,3)-xs(2,3)*xs(3,1))
140     &   +xs(1,3)*(xs(2,1)*xs(3,2)-xs(2,2)*xs(3,1))
141!
142      if(iflag.eq.2) return
143!
144!     computation of the global derivative of the local coordinates
145!     (xsi) (inversion of xs)
146!
147      xsi(1,1)=(xs(2,2)*xs(3,3)-xs(3,2)*xs(2,3))/xsj
148      xsi(1,2)=(xs(1,3)*xs(3,2)-xs(1,2)*xs(3,3))/xsj
149      xsi(1,3)=(xs(1,2)*xs(2,3)-xs(2,2)*xs(1,3))/xsj
150      xsi(2,1)=(xs(2,3)*xs(3,1)-xs(2,1)*xs(3,3))/xsj
151      xsi(2,2)=(xs(1,1)*xs(3,3)-xs(3,1)*xs(1,3))/xsj
152      xsi(2,3)=(xs(1,3)*xs(2,1)-xs(1,1)*xs(2,3))/xsj
153      xsi(3,1)=(xs(2,1)*xs(3,2)-xs(3,1)*xs(2,2))/xsj
154      xsi(3,2)=(xs(1,2)*xs(3,1)-xs(1,1)*xs(3,2))/xsj
155      xsi(3,3)=(xs(1,1)*xs(2,2)-xs(2,1)*xs(1,2))/xsj
156!
157!     computation of the global derivatives of the shape functions
158!
159      do k=1,15
160        do j=1,3
161          sh(j)=shp(1,k)*xsi(1,j)+shp(2,k)*xsi(2,j)+shp(3,k)*xsi(3,j)
162        enddo
163        do j=1,3
164          shp(j,k)=sh(j)
165        enddo
166      enddo
167!
168      return
169      end
170