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 dual shape funciton \f$ shape(\xi,\eta) \f$
21!
22!  [in] xi		xi-coordinate
23!  [in] et		eta-coordinate
24!  [in] xl		local node coordinates
25!  [out] xsj		jacobian vector
26!  [out] xs		local derivative of the global coordinates
27!  [out] shp		evaluated shape functions and derivatives
28!  [in] ns		current slave face
29!  [in] pslavdual 	(:,i)coefficients \f$ \alpha_{ij}\f$,
30!                       \f$ 1,j=1,..8\f$ for dual shape functions for face i
31!  [in] iflag		flag indicating what to compute
32!
33      subroutine dualshape4q(xi,et,xl,xsj,xs,shp,ns,pslavdual,iflag)
34!
35!     iflag=2: calculate the value of the shape functions,
36!              their derivatives w.r.t. the local coordinates
37!              and the Jacobian vector (local normal to the
38!              surface)
39!     iflag=3: calculate the value of the shape functions, the
40!              value of their derivatives w.r.t. the global
41!              coordinates and the Jacobian vector (local normal
42!              to the surface)
43!
44      implicit none
45!
46      integer i,j,k,iflag,ns
47!
48      real*8 shp(7,8),xs(3,7),xsi(2,3),xl(3,8),sh(3),xsj(3)
49!
50      real*8 xi,et,pslavdual(64,*)
51!
52!
53!
54!     shape functions and their glocal derivatives for an element
55!     described with two local parameters and three global ones.
56!
57!     Caution: derivatives and exspecially jacobian for standard
58!     basis functions are given
59!     needed for consistent integration
60!
61!     local derivatives of the shape functions: xi-derivative
62!
63      shp(1,1)=-(1.d0-et)/4.d0
64      shp(1,2)=(1.d0-et)/4.d0
65      shp(1,3)=(1.d0+et)/4.d0
66      shp(1,4)=-(1.d0+et)/4.d0
67!
68!     local derivatives of the shape functions: eta-derivative
69!
70      shp(2,1)=-(1.d0-xi)/4.d0
71      shp(2,2)=-(1.d0+xi)/4.d0
72      shp(2,3)=(1.d0+xi)/4.d0
73      shp(2,4)=(1.d0-xi)/4.d0
74!
75!     standard shape functions
76!
77      shp(3,1)=(1.d0-xi)*(1.d0-et)/4.d0
78      shp(3,2)=(1.d0+xi)*(1.d0-et)/4.d0
79      shp(3,3)=(1.d0+xi)*(1.d0+et)/4.d0
80      shp(3,4)=(1.d0-xi)*(1.d0+et)/4.d0
81!
82!     Dual shape functions
83!     with Mass Matrix pslavdual
84!
85      do i=1,4
86         shp(4,i)=0.0
87         do j=1,4
88            shp(4,i)=shp(4,i)+pslavdual((i-1)*8+j,ns)*shp(3,j)
89         enddo
90      enddo
91!
92!     computation of the local derivative of the global coordinates
93!     (xs)
94!
95      do i=1,3
96         do j=1,2
97            xs(i,j)=0.d0
98            do k=1,4
99               xs(i,j)=xs(i,j)+xl(i,k)*shp(j,k)
100            enddo
101         enddo
102      enddo
103!
104!     computation of the jacobian vector
105!
106      xsj(1)=xs(2,1)*xs(3,2)-xs(3,1)*xs(2,2)
107      xsj(2)=xs(1,2)*xs(3,1)-xs(3,2)*xs(1,1)
108      xsj(3)=xs(1,1)*xs(2,2)-xs(2,1)*xs(1,2)
109!
110      if(iflag.eq.2) return
111!
112!     computation of the global derivative of the local coordinates
113!     (xsi) (inversion of xs)
114!
115      xsi(1,1)=xs(2,2)/xsj(3)
116      xsi(2,1)=-xs(2,1)/xsj(3)
117      xsi(1,2)=-xs(1,2)/xsj(3)
118      xsi(2,2)=xs(1,1)/xsj(3)
119      xsi(1,3)=-xs(2,2)/xsj(1)
120      xsi(2,3)=xs(2,1)/xsj(1)
121!
122!     computation of the global derivatives of the shape functions
123!
124      do k=1,4
125         do j=1,3
126            sh(j)=shp(1,k)*xsi(1,j)+shp(2,k)*xsi(2,j)
127         enddo
128         do j=1,3
129            shp(j,k)=sh(j)
130         enddo
131      enddo
132!
133      return
134      end
135
136