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!
20!     function to evaluate transformed shape funciton \f$ shape(\xi,\eta) \f$
21!     for quad-quad mortar method, see phd-thesis Sitzmann Chapter 4.1.
22!
23!     Author: Saskia Sitzmann
24!
25!  [in] xi		xi-coordinate
26!  [in] et		eta-coordinate
27!  [in] xl		local node coordinates
28!  [out] xsj		jacobian vector
29!  [out] xs		local derivative of the global coordinates
30!  [out] shp		evaluated shape functions and derivatives
31!  [in] iflag		flag indicating what to compute
32!
33      subroutine shape8qtilde(xi,et,xl,xsj,xs,shp,iflag)
34!
35!     shape functions and derivatives for a 8-node quadratic
36!     isoparametric quadrilateral element. -1<=xi,et<=1
37!
38!     iflag=1: calculate only the value of the shape functions
39!     iflag=2: calculate the value of the shape functions,
40!              their derivatives w.r.t. the local coordinates
41!              and the Jacobian vector (local normal to the
42!              surface)
43!     iflag=3: calculate the value of the shape functions, the
44!              value of their derivatives w.r.t. the global
45!              coordinates and the Jacobian vector (local normal
46!              to the surface)
47!     iflag=4: calculate the value of the shape functions, the
48!              value of their 1st and 2nd order derivatives
49!              w.r.t. the local coordinates, the Jacobian vector
50!              (local normal to the surface)
51!
52      implicit none
53!
54      integer i,j,k,iflag
55!
56      real*8 shp(7,8),xs(3,7),xsi(2,3),xl(3,8),sh(3),xsj(3),xi,
57     &     et,shpold(7,8),alpha
58!
59!
60!
61!     shape functions and their glocal derivatives for an element
62!     described with two local parameters and three global ones.
63!
64      alpha=1.0/5.0
65!
66!     shape functions
67!
68      shpold(4,1)=(1.d0-xi)*(1.d0-et)*(-xi-et-1.d0)/4.d0
69      shpold(4,2)=(1.d0+xi)*(1.d0-et)*(xi-et-1.d0)/4.d0
70      shpold(4,3)=(1.d0+xi)*(1.d0+et)*(xi+et-1.d0)/4.d0
71      shpold(4,4)=(1.d0-xi)*(1.d0+et)*(-xi+et-1.d0)/4.d0
72      shpold(4,5)=(1.d0-xi*xi)*(1.d0-et)/2.d0
73      shpold(4,6)=(1.d0+xi)*(1.d0-et*et)/2.d0
74      shpold(4,7)=(1.d0-xi*xi)*(1.d0+et)/2.d0
75      shpold(4,8)=(1.d0-xi)*(1.d0-et*et)/2.d0
76!
77      shp(4,1)=1.0*shpold(4,1)+alpha*shpold(4,5)+alpha*shpold(4,8)
78      shp(4,2)=1.0*shpold(4,2)+alpha*shpold(4,5)+alpha*shpold(4,6)
79      shp(4,3)=1.0*shpold(4,3)+alpha*shpold(4,6)+alpha*shpold(4,7)
80      shp(4,4)=1.0*shpold(4,4)+alpha*shpold(4,7)+alpha*shpold(4,8)
81      shp(4,5)=(1.0-2.0*alpha)*shpold(4,5)
82      shp(4,6)=(1.0-2.0*alpha)*shpold(4,6)
83      shp(4,7)=(1.0-2.0*alpha)*shpold(4,7)
84      shp(4,8)=(1.0-2.0*alpha)*shpold(4,8)
85!
86!     Caution: derivatives and exspecially jacobian for untransformed
87!     basis functions are given
88!     needed for consistent integration
89!
90      if(iflag.eq.1) return
91!
92!     local derivatives of the shape functions: xi-derivative
93!
94      shp(1,1)=(1.d0-et)*(2.d0*xi+et)/4.d0
95      shp(1,2)=(1.d0-et)*(2.d0*xi-et)/4.d0
96      shp(1,3)=(1.d0+et)*(2.d0*xi+et)/4.d0
97      shp(1,4)=(1.d0+et)*(2.d0*xi-et)/4.d0
98      shp(1,5)=-xi*(1.d0-et)
99      shp(1,6)=(1.d0-et*et)/2.d0
100      shp(1,7)=-xi*(1.d0+et)
101      shp(1,8)=-(1.d0-et*et)/2.d0
102!
103!     local derivatives of the shape functions: eta-derivative
104!
105      shp(2,1)=(1.d0-xi)*(2.d0*et+xi)/4.d0
106      shp(2,2)=(1.d0+xi)*(2.d0*et-xi)/4.d0
107      shp(2,3)=(1.d0+xi)*(2.d0*et+xi)/4.d0
108      shp(2,4)=(1.d0-xi)*(2.d0*et-xi)/4.d0
109      shp(2,5)=-(1.d0-xi*xi)/2.d0
110      shp(2,6)=-et*(1.d0+xi)
111      shp(2,7)=(1.d0-xi*xi)/2.d0
112      shp(2,8)=-et*(1.d0-xi)
113!
114!     computation of the local derivative of the global coordinates
115!     (xs)
116!
117      do i=1,3
118         do j=1,2
119            xs(i,j)=0.d0
120            do k=1,8
121               xs(i,j)=xs(i,j)+xl(i,k)*shp(j,k)
122            enddo
123         enddo
124      enddo
125!
126!     computation of the jacobian vector
127!
128      xsj(1)=xs(2,1)*xs(3,2)-xs(3,1)*xs(2,2)
129      xsj(2)=xs(1,2)*xs(3,1)-xs(3,2)*xs(1,1)
130      xsj(3)=xs(1,1)*xs(2,2)-xs(2,1)*xs(1,2)
131!
132      if(iflag.eq.3) then
133!
134!     computation of the global derivative of the local coordinates
135!     (xsi) (inversion of xs)
136!
137         if(dabs(xsj(3)).gt.1.d-10) then
138            xsi(1,1)=xs(2,2)/xsj(3)
139            xsi(2,2)=xs(1,1)/xsj(3)
140            xsi(1,2)=-xs(1,2)/xsj(3)
141            xsi(2,1)=-xs(2,1)/xsj(3)
142            if(dabs(xsj(2)).gt.1.d-10) then
143               xsi(2,3)=xs(1,1)/(-xsj(2))
144               xsi(1,3)=-xs(1,2)/(-xsj(2))
145            elseif(dabs(xsj(1)).gt.1.d-10) then
146               xsi(2,3)=xs(2,1)/xsj(1)
147               xsi(1,3)=-xs(2,2)/xsj(1)
148            else
149               xsi(2,3)=0.d0
150               xsi(1,3)=0.d0
151            endif
152         elseif(dabs(xsj(2)).gt.1.d-10) then
153            xsi(1,1)=xs(3,2)/(-xsj(2))
154            xsi(2,3)=xs(1,1)/(-xsj(2))
155            xsi(1,3)=-xs(1,2)/(-xsj(2))
156            xsi(2,1)=-xs(3,1)/(-xsj(2))
157            if(dabs(xsj(1)).gt.1.d-10) then
158               xsi(1,2)=xs(3,2)/xsj(1)
159               xsi(2,2)=-xs(3,1)/xsj(1)
160            else
161               xsi(1,2)=0.d0
162               xsi(2,2)=0.d0
163            endif
164         else
165            xsi(1,2)=xs(3,2)/xsj(1)
166            xsi(2,3)=xs(2,1)/xsj(1)
167            xsi(1,3)=-xs(2,2)/xsj(1)
168            xsi(2,2)=-xs(3,1)/xsj(1)
169            xsi(1,1)=0.d0
170            xsi(2,1)=0.d0
171         endif
172!
173!     computation of the global derivatives of the shape functions
174!
175         do k=1,8
176            do j=1,3
177               sh(j)=shp(1,k)*xsi(1,j)+shp(2,k)*xsi(2,j)
178            enddo
179            do j=1,3
180               shp(j,k)=sh(j)
181            enddo
182         enddo
183!
184      elseif(iflag.eq.4) then
185!
186!     local 2nd order derivatives of the shape functions: xi,xi-derivative
187!
188         shp(5,1)=(1.d0-et)/2.d0
189         shp(5,2)=(1.d0-et)/2.d0
190         shp(5,3)=(1.d0+et)/2.d0
191         shp(5,4)=(1.d0+et)/2.d0
192         shp(5,5)=-(1.d0-et)
193         shp(5,6)=0.d0
194         shp(5,7)=-(1.d0+et)
195         shp(5,8)=0.d0
196!
197         shp(5,1)=1.0*shpold(5,1)+alpha*shpold(5,5)+alpha*shpold(5,8)
198         shp(5,2)=1.0*shpold(5,2)+alpha*shpold(5,5)+alpha*shpold(5,6)
199         shp(5,3)=1.0*shpold(5,3)+alpha*shpold(5,6)+alpha*shpold(5,7)
200         shp(5,4)=1.0*shpold(5,4)+alpha*shpold(5,7)+alpha*shpold(5,8)
201         shp(5,5)=(1.0-2.0*alpha)*shpold(5,5)
202         shp(5,6)=(1.0-2.0*alpha)*shpold(5,6)
203         shp(5,7)=(1.0-2.0*alpha)*shpold(5,7)
204         shp(5,8)=(1.0-2.0*alpha)*shpold(5,8)
205!
206!     local 2nd order derivatives of the shape functions: xi,eta-derivative
207!
208         shp(6,1)=(1.d0-2.d0*(xi+et))/4.d0
209         shp(6,2)=(-1.d0-2.d0*(xi-et))/4.d0
210         shp(6,3)=(1.d0+2.d0*(xi+et))/4.d0
211         shp(6,4)=(-1.d0-2.d0*(xi+et))/4.d0
212         shp(6,5)=xi
213         shp(6,6)=-et
214         shp(6,7)=-xi
215         shp(6,8)=et
216!
217         shp(6,1)=1.0*shpold(6,1)+alpha*shpold(6,5)+alpha*shpold(6,8)
218         shp(6,2)=1.0*shpold(6,2)+alpha*shpold(6,5)+alpha*shpold(6,6)
219         shp(6,3)=1.0*shpold(6,3)+alpha*shpold(6,6)+alpha*shpold(6,7)
220         shp(6,4)=1.0*shpold(6,4)+alpha*shpold(6,7)+alpha*shpold(6,8)
221         shp(6,5)=(1.0-2.0*alpha)*shpold(6,5)
222         shp(6,6)=(1.0-2.0*alpha)*shpold(6,6)
223         shp(6,7)=(1.0-2.0*alpha)*shpold(6,7)
224         shp(6,8)=(1.0-2.0*alpha)*shpold(6,8)
225!
226!     local 2nd order derivatives of the shape functions: eta,eta-derivative
227!
228         shp(7,1)=(1.d0-xi)/2.d0
229         shp(7,2)=(1.d0+xi)/2.d0
230         shp(7,3)=(1.d0+xi)/2.d0
231         shp(7,4)=(1.d0-xi)/2.d0
232         shp(7,5)=0.d0
233         shp(7,6)=-(1.d0+xi)
234         shp(7,7)=0.d0
235         shp(7,8)=-(1.d0-xi)
236!
237         shp(7,1)=1.0*shpold(7,1)+alpha*shpold(7,5)+alpha*shpold(7,8)
238         shp(7,2)=1.0*shpold(7,2)+alpha*shpold(7,5)+alpha*shpold(7,6)
239         shp(7,3)=1.0*shpold(7,3)+alpha*shpold(7,6)+alpha*shpold(7,7)
240         shp(7,4)=1.0*shpold(7,4)+alpha*shpold(7,7)+alpha*shpold(7,8)
241         shp(7,5)=(1.0-2.0*alpha)*shpold(7,5)
242         shp(7,6)=(1.0-2.0*alpha)*shpold(7,6)
243         shp(7,7)=(1.0-2.0*alpha)*shpold(7,7)
244         shp(7,8)=(1.0-2.0*alpha)*shpold(7,8)
245!
246!     computation of the local 2nd derivatives of the global coordinates
247!     (xs)
248!
249         do i=1,3
250            do j=5,7
251               xs(i,j)=0.d0
252               do k=1,8
253                  xs(i,j)=xs(i,j)+xl(i,k)*shp(j,k)
254               enddo
255            enddo
256         enddo
257      endif
258!
259      return
260      end
261
262