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 shape funciton dual shape funciton for
21!     quad-quad mortar method \f$ shape(\xi,\eta) \f$
22!     see phd-thesis Sitzmann Chapter 4.1.
23!
24!     Author: Saskia Sitzmann
25!
26!  [in] xi		xi-coordinate
27!  [in] et		eta-coordinate
28!  [in] xl		local node coordinates
29!  [out] xsj		jacobian vector
30!  [out] xs		local derivative of the global coordinates
31!  [out] shp		evaluated shape functions and derivatives
32!  [in] ns		current slave face
33!  [in] pslavdual 	(:,i)coefficients \f$ \alpha_{ij}\f$,
34!                       \f$ 1,j=1,..8\f$ for dual shape functions for face i
35!  [in] iflag		flag indicating what to compute
36!
37      subroutine dualshape8qtilde(xi,et,xl,xsj,xs,shp,ns,
38     &     pslavdual,iflag)
39!
40!     shape functions and derivatives for a 8-node quadratic
41!     isoparametric quadrilateral element. -1<=xi,et<=1
42!
43!     iflag=2: calculate the value of the shape functions,
44!     their derivatives w.r.t. the local coordinates
45!              and the Jacobian vector (local normal to the
46!              surface)
47!     iflag=3: calculate the value of the shape functions, the
48!              value of their derivatives w.r.t. the global
49!              coordinates and the Jacobian vector (local normal
50!              to the surface)
51!
52      implicit none
53!
54      integer i,j,k,iflag,ns
55!
56      real*8 shp(7,8),xs(3,7),xsi(2,3),xl(3,8),sh(3),xsj(3),
57     &     shpold(7,8),alpha,pslavdual(64,*)
58!
59      real*8 xi,et
60!
61!
62!
63!     shape functions and their glocal derivatives for an element
64!     described with two local parameters and three global ones.
65!
66      alpha=1.0/5.0
67!
68!     standard shape functions
69!
70!     shape functions
71!
72      shpold(4,1)=(1.d0-xi)*(1.d0-et)*(-xi-et-1.d0)/4.d0
73      shpold(4,2)=(1.d0+xi)*(1.d0-et)*(xi-et-1.d0)/4.d0
74      shpold(4,3)=(1.d0+xi)*(1.d0+et)*(xi+et-1.d0)/4.d0
75      shpold(4,4)=(1.d0-xi)*(1.d0+et)*(-xi+et-1.d0)/4.d0
76      shpold(4,5)=(1.d0-xi*xi)*(1.d0-et)/2.d0
77      shpold(4,6)=(1.d0+xi)*(1.d0-et*et)/2.d0
78      shpold(4,7)=(1.d0-xi*xi)*(1.d0+et)/2.d0
79      shpold(4,8)=(1.d0-xi)*(1.d0-et*et)/2.d0
80!
81!     transformed standard shape functions
82!
83      shp(3,1)=1.0*shpold(4,1)+alpha*shpold(4,5)+alpha*shpold(4,8)
84      shp(3,2)=1.0*shpold(4,2)+alpha*shpold(4,5)+alpha*shpold(4,6)
85      shp(3,3)=1.0*shpold(4,3)+alpha*shpold(4,6)+alpha*shpold(4,7)
86      shp(3,4)=1.0*shpold(4,4)+alpha*shpold(4,7)+alpha*shpold(4,8)
87      shp(3,5)=(1.0-2.0*alpha)*shpold(4,5)
88      shp(3,6)=(1.0-2.0*alpha)*shpold(4,6)
89      shp(3,7)=(1.0-2.0*alpha)*shpold(4,7)
90      shp(3,8)=(1.0-2.0*alpha)*shpold(4,8)
91!
92!     Caution: derivatives and exspecially jacobian for untransformed
93!     basis functions are given
94!     needed for consistent integration
95!
96!     local derivatives of the shape functions: xi-derivative
97!
98      shp(1,1)=(1.d0-et)*(2.d0*xi+et)/4.d0
99      shp(1,2)=(1.d0-et)*(2.d0*xi-et)/4.d0
100      shp(1,3)=(1.d0+et)*(2.d0*xi+et)/4.d0
101      shp(1,4)=(1.d0+et)*(2.d0*xi-et)/4.d0
102      shp(1,5)=-xi*(1.d0-et)
103      shp(1,6)=(1.d0-et*et)/2.d0
104      shp(1,7)=-xi*(1.d0+et)
105      shp(1,8)=-(1.d0-et*et)/2.d0
106!
107!     local derivatives of the shape functions: eta-derivative
108!
109      shp(2,1)=(1.d0-xi)*(2.d0*et+xi)/4.d0
110      shp(2,2)=(1.d0+xi)*(2.d0*et-xi)/4.d0
111      shp(2,3)=(1.d0+xi)*(2.d0*et+xi)/4.d0
112      shp(2,4)=(1.d0-xi)*(2.d0*et-xi)/4.d0
113      shp(2,5)=-(1.d0-xi*xi)/2.d0
114      shp(2,6)=-et*(1.d0+xi)
115      shp(2,7)=(1.d0-xi*xi)/2.d0
116      shp(2,8)=-et*(1.d0-xi)
117!
118!     Dual shape functions with transformed std basis
119!
120      do i=1,8
121         shp(4,i)=0.0
122         do j=1,8
123            shp(4,i)=shp(4,i)+pslavdual((i-1)*8+j,ns)*shp(3,j)
124         enddo
125      enddo
126!
127!     computation of the local derivative of the global coordinates
128!     (xs)
129!
130      do i=1,3
131         do j=1,2
132            xs(i,j)=0.d0
133            do k=1,8
134               xs(i,j)=xs(i,j)+xl(i,k)*shp(j,k)
135            enddo
136         enddo
137      enddo
138!
139!     computation of the jacobian vector
140!
141      xsj(1)=xs(2,1)*xs(3,2)-xs(3,1)*xs(2,2)
142      xsj(2)=xs(1,2)*xs(3,1)-xs(3,2)*xs(1,1)
143      xsj(3)=xs(1,1)*xs(2,2)-xs(2,1)*xs(1,2)
144!
145      if(iflag.eq.2) return
146!
147!     computation of the global derivative of the local coordinates
148!     (xsi) (inversion of xs)
149!
150      if(dabs(xsj(3)).gt.1.d-10) then
151         xsi(1,1)=xs(2,2)/xsj(3)
152         xsi(2,2)=xs(1,1)/xsj(3)
153         xsi(1,2)=-xs(1,2)/xsj(3)
154         xsi(2,1)=-xs(2,1)/xsj(3)
155         if(dabs(xsj(2)).gt.1.d-10) then
156            xsi(2,3)=xs(1,1)/xsj(2)
157            xsi(1,3)=-xs(1,2)/xsj(2)
158         elseif(dabs(xsj(1)).gt.1.d-10) then
159            xsi(2,3)=xs(2,1)/xsj(1)
160            xsi(1,3)=-xs(2,2)/xsj(1)
161         else
162            xsi(2,3)=0.d0
163            xsi(1,3)=0.d0
164         endif
165      elseif(dabs(xsj(2)).gt.1.d-10) then
166         xsi(1,1)=xs(3,2)/xsj(2)
167         xsi(2,3)=xs(1,1)/xsj(2)
168         xsi(1,3)=-xs(1,2)/xsj(2)
169         xsi(2,1)=-xs(3,1)/xsj(2)
170         if(dabs(xsj(1)).gt.1.d-10) then
171            xsi(1,2)=xs(3,2)/xsj(1)
172            xsi(2,2)=-xs(3,1)/xsj(1)
173         else
174            xsi(1,2)=0.d0
175            xsi(2,2)=0.d0
176         endif
177      else
178         xsi(1,2)=xs(3,2)/xsj(1)
179         xsi(2,3)=xs(2,1)/xsj(1)
180         xsi(1,3)=-xs(2,2)/xsj(1)
181         xsi(2,2)=-xs(3,1)/xsj(1)
182         xsi(1,1)=0.d0
183         xsi(2,1)=0.d0
184      endif
185!
186!     computation of the global derivatives of the shape functions
187!
188      do k=1,8
189         do j=1,3
190            sh(j)=shp(1,k)*xsi(1,j)+shp(2,k)*xsi(2,j)
191         enddo
192         do j=1,3
193            shp(j,k)=sh(j)
194         enddo
195      enddo
196!
197      return
198      end
199
200