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 shape9q(xi,et,xl,xsj,xs,shp,iflag)
20!
21!     shape functions and derivatives for a 9-node quadratic
22!     isoparametric quadrilateral element. -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,
26!              their derivatives w.r.t. the local coordinates
27!              and the Jacobian vector (local normal to the
28!              surface)
29!     iflag=3: calculate the value of the shape functions, the
30!              value of their derivatives w.r.t. the global
31!              coordinates and the Jacobian vector (local normal
32!              to the surface)
33!     iflag=4: calculate the value of the shape functions, the
34!              value of their 1st and 2nd order derivatives
35!              w.r.t. the local coordinates, the Jacobian vector
36!              (local normal to the surface)
37!
38      implicit none
39!
40      integer i,j,k,iflag
41!
42      real*8 shp(7,9),xs(3,7),xsi(2,3),xl(3,9),sh(3),xsj(3),xi,et,
43     &  fxi1,fxi2,fxi3,fet1,fet2,fet3,dfxi1,dfxi2,dfxi3,dfet1,dfet2,
44     &  dfet3,ddfxi1,ddfxi2,ddfxi3,ddfet1,ddfet2,ddfet3
45!
46!
47!
48!     shape functions and their glocal derivatives for an element
49!     described with two local parameters and three global ones.
50!
51!     shape functions in one dimension
52!
53      fxi1=xi*(xi-1.d0)/2.d0
54      fxi2=(1.d0-xi)*(1.d0+xi)
55      fxi3=xi*(xi+1.d0)/2.d0
56!
57      fet1=et*(et-1.d0)/2.d0
58      fet2=(1.d0-et)*(1.d0+et)
59      fet3=et*(et+1.d0)/2.d0
60!
61!     shape functions
62!
63      shp(4,1)=fxi1*fet1
64      shp(4,2)=fxi3*fet1
65      shp(4,3)=fxi3*fet3
66      shp(4,4)=fxi1*fet3
67      shp(4,5)=fxi2*fet1
68      shp(4,6)=fxi3*fet2
69      shp(4,7)=fxi2*fet3
70      shp(4,8)=fxi1*fet2
71      shp(4,9)=fxi2*fet2
72!
73      if(iflag.eq.1) return
74!
75!     derivative of the shape functions in one dimension
76!
77      dfxi1=(2.d0*xi-1.d0)/2.d0
78      dfxi2=-2.d0*xi
79      dfxi3=(2.d0*xi+1.d0)/2.d0
80!
81      dfet1=(2.d0*et-1.d0)/2.d0
82      dfet2=-2.d0*et
83      dfet3=(2.d0*et+1.d0)/2.d0
84!
85!     local derivatives of the shape functions: xi-derivative
86!
87      shp(1,1)=dfxi1*fet1
88      shp(1,2)=dfxi3*fet1
89      shp(1,3)=dfxi3*fet3
90      shp(1,4)=dfxi1*fet3
91      shp(1,5)=dfxi2*fet1
92      shp(1,6)=dfxi3*fet2
93      shp(1,7)=dfxi2*fet3
94      shp(1,8)=dfxi1*fet2
95      shp(1,9)=dfxi2*fet2
96!
97!     local derivatives of the shape functions: eta-derivative
98!
99      shp(2,1)=fxi1*dfet1
100      shp(2,2)=fxi3*dfet1
101      shp(2,3)=fxi3*dfet3
102      shp(2,4)=fxi1*dfet3
103      shp(2,5)=fxi2*dfet1
104      shp(2,6)=fxi3*dfet2
105      shp(2,7)=fxi2*dfet3
106      shp(2,8)=fxi1*dfet2
107      shp(2,9)=fxi2*dfet2
108!
109!     computation of the local derivative of the global coordinates
110!     (xs)
111!
112      do i=1,3
113        do j=1,2
114          xs(i,j)=0.d0
115          do k=1,9
116            xs(i,j)=xs(i,j)+xl(i,k)*shp(j,k)
117          enddo
118        enddo
119      enddo
120!
121!     computation of the jacobian vector
122!
123      xsj(1)=xs(2,1)*xs(3,2)-xs(3,1)*xs(2,2)
124      xsj(2)=xs(1,2)*xs(3,1)-xs(3,2)*xs(1,1)
125      xsj(3)=xs(1,1)*xs(2,2)-xs(2,1)*xs(1,2)
126!
127      if(iflag.eq.3) then
128!
129!     computation of the global derivative of the local coordinates
130!     (xsi) (inversion of xs)
131!
132         if(dabs(xsj(3)).gt.1.d-10) then
133            xsi(1,1)=xs(2,2)/xsj(3)
134            xsi(2,2)=xs(1,1)/xsj(3)
135            xsi(1,2)=-xs(1,2)/xsj(3)
136            xsi(2,1)=-xs(2,1)/xsj(3)
137            if(dabs(xsj(2)).gt.1.d-10) then
138               xsi(2,3)=xs(1,1)/(-xsj(2))
139               xsi(1,3)=-xs(1,2)/(-xsj(2))
140            elseif(dabs(xsj(1)).gt.1.d-10) then
141               xsi(2,3)=xs(2,1)/xsj(1)
142               xsi(1,3)=-xs(2,2)/xsj(1)
143            else
144               xsi(2,3)=0.d0
145               xsi(1,3)=0.d0
146            endif
147         elseif(dabs(xsj(2)).gt.1.d-10) then
148            xsi(1,1)=xs(3,2)/(-xsj(2))
149            xsi(2,3)=xs(1,1)/(-xsj(2))
150            xsi(1,3)=-xs(1,2)/(-xsj(2))
151            xsi(2,1)=-xs(3,1)/(-xsj(2))
152            if(dabs(xsj(1)).gt.1.d-10) then
153               xsi(1,2)=xs(3,2)/xsj(1)
154               xsi(2,2)=-xs(3,1)/xsj(1)
155            else
156               xsi(1,2)=0.d0
157               xsi(2,2)=0.d0
158            endif
159         else
160            xsi(1,2)=xs(3,2)/xsj(1)
161            xsi(2,3)=xs(2,1)/xsj(1)
162            xsi(1,3)=-xs(2,2)/xsj(1)
163            xsi(2,2)=-xs(3,1)/xsj(1)
164            xsi(1,1)=0.d0
165            xsi(2,1)=0.d0
166         endif
167!
168!     computation of the global derivatives of the shape functions
169!
170         do k=1,9
171            do j=1,3
172               sh(j)=shp(1,k)*xsi(1,j)+shp(2,k)*xsi(2,j)
173            enddo
174            do j=1,3
175               shp(j,k)=sh(j)
176            enddo
177         enddo
178!
179      elseif(iflag.eq.4) then
180!
181!     second derivative of the shape functions in one dimension
182!
183         ddfxi1=1.d0
184         ddfxi2=-2.d0
185         ddfxi3=1.d0
186!
187         ddfet1=1.d0
188         ddfet2=-2.d0
189         ddfet3=1.d0
190!
191!     local 2nd order derivatives of the shape functions: xi,xi-derivative
192!
193         shp(5,1)=ddfxi1*fet1
194         shp(5,2)=ddfxi3*fet1
195         shp(5,3)=ddfxi3*fet3
196         shp(5,4)=ddfxi1*fet3
197         shp(5,5)=ddfxi2*fet1
198         shp(5,6)=ddfxi3*fet2
199         shp(5,7)=ddfxi2*fet3
200         shp(5,8)=ddfxi1*fet2
201         shp(5,9)=ddfxi2*fet2
202!
203!     local 2nd order derivatives of the shape functions: xi,eta-derivative
204!
205         shp(6,1)=dfxi1*dfet1
206         shp(6,2)=dfxi3*dfet1
207         shp(6,3)=dfxi3*dfet3
208         shp(6,4)=dfxi1*dfet3
209         shp(6,5)=dfxi2*dfet1
210         shp(6,6)=dfxi3*dfet2
211         shp(6,7)=dfxi2*dfet3
212         shp(6,8)=dfxi1*dfet2
213         shp(6,9)=dfxi2*dfet2
214!
215!     local 2nd order derivatives of the shape functions: eta,eta-derivative
216!
217         shp(7,1)=fxi1*ddfet1
218         shp(7,2)=fxi3*ddfet1
219         shp(7,3)=fxi3*ddfet3
220         shp(7,4)=fxi1*ddfet3
221         shp(7,5)=fxi2*ddfet1
222         shp(7,6)=fxi3*ddfet2
223         shp(7,7)=fxi2*ddfet3
224         shp(7,8)=fxi1*ddfet2
225         shp(7,9)=fxi2*ddfet2
226!
227!     computation of the local 2nd derivatives of the global coordinates
228!     (xs)
229!
230         do i=1,3
231            do j=5,7
232               xs(i,j)=0.d0
233               do k=1,9
234                  xs(i,j)=xs(i,j)+xl(i,k)*shp(j,k)
235               enddo
236            enddo
237         enddo
238      endif
239!
240      return
241      end
242