1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
5!> \brief  This module contains functions for interpolation in 9 node
6!!   quadrilateral element
7module shape_quad9n
8  !####################################################################
9
10  implicit none
11
12  !--------------------------------------------------------------------
13
14  integer, parameter, private :: kreal = kind( 0.0D0 )
15
16  !--------------------------------------------------------------------
17
18contains
19
20
21  !####################################################################
22  subroutine ShapeFunc_quad9n(lcoord, func)
23    !####################################################################
24
25    real(kind = kreal), intent(in)  :: lcoord(2)
26    real(kind = kreal), intent(out) :: func(9)
27
28    !--------------------------------------------------------------------
29
30    real(kind = kreal) :: xi(9), eta(9)
31    real(kind = kreal) :: n_xi(9), n_eta(9)
32
33    integer :: na
34
35    !--------------------------------------------------------------------
36
37    ! xi-coordinate at a node in a local element
38    xi(1)  = -1.0D0
39    xi(2)  =  1.0D0
40    xi(3)  =  1.0D0
41    xi(4)  = -1.0D0
42    xi(5)  =  0.0D0
43    xi(6)  =  1.0D0
44    xi(7)  =  0.0D0
45    xi(8)  = -1.0D0
46    xi(9)  =  0.0D0
47    ! eta-coordinate at a node in a local element
48    eta(1) = -1.0D0
49    eta(2) = -1.0D0
50    eta(3) =  1.0D0
51    eta(4) =  1.0D0
52    eta(5) = -1.0D0
53    eta(6) =  0.0D0
54    eta(7) =  1.0D0
55    eta(8) =  0.0D0
56    eta(9) =  0.0D0
57
58    !--------------------------------------------------------------------
59
60    do na = 1, 9
61
62      n_xi(na)  = ( 0.5D0*xi(na) *lcoord(1) )    &
63        *( 1.0D0+xi(na) *lcoord(1) )   &
64        +( 1.0D0-xi(na) *xi(na)  )      &
65        *( 1.0D0-lcoord(1)*lcoord(1) )
66      n_eta(na) = ( 0.5D0*eta(na)*lcoord(2) )    &
67        *( 1.0D0+eta(na)*lcoord(2) )   &
68        +( 1.0D0-eta(na)*eta(na) )      &
69        *( 1.0D0-lcoord(2)*lcoord(2) )
70
71    end do
72
73    !--------------------------------------------------------------------
74
75    do na = 1, 9
76
77      func(na) = n_xi(na)*n_eta(na)
78
79    end do
80
81    !--------------------------------------------------------------------
82
83    return
84
85    !####################################################################
86  end subroutine ShapeFunc_quad9n
87  !####################################################################
88
89  !####################################################################
90  subroutine ShapeDeriv_quad9n(lcoord, func)
91    !####################################################################
92
93    real(kind = kreal), intent(in)  :: lcoord(2)
94    real(kind = kreal), intent(out) :: func(9, 2)
95
96    !--------------------------------------------------------------------
97
98    real(kind = kreal) :: xi(9), eta(9)
99    real(kind = kreal) :: n_xi(9), n_eta(9)
100    real(kind = kreal) :: dn_xi(9), dn_eta(9)
101
102    integer :: na
103
104    !--------------------------------------------------------------------
105
106    ! xi-coordinate at a node in a local element
107    xi(1)  = -1.0D0
108    xi(2)  =  1.0D0
109    xi(3)  =  1.0D0
110    xi(4)  = -1.0D0
111    xi(5)  =  0.0D0
112    xi(6)  =  1.0D0
113    xi(7)  =  0.0D0
114    xi(8)  = -1.0D0
115    xi(9)  =  0.0D0
116    ! eta-coordinate at a node in a local element
117    eta(1) = -1.0D0
118    eta(2) = -1.0D0
119    eta(3) =  1.0D0
120    eta(4) =  1.0D0
121    eta(5) = -1.0D0
122    eta(6) =  0.0D0
123    eta(7) =  1.0D0
124    eta(8) =  0.0D0
125    eta(9) =  0.0D0
126
127    !--------------------------------------------------------------------
128
129    do na = 1, 9
130
131      n_xi(na)  = ( 0.5D0*xi(na) *lcoord(1) )    &
132        *( 1.0D0+xi(na) *lcoord(1) )   &
133        +( 1.0D0-xi(na) *xi(na)  )      &
134        *( 1.0D0-lcoord(1)*lcoord(1) )
135      n_eta(na) = ( 0.5D0*eta(na)*lcoord(2) )    &
136        *( 1.0D0+eta(na)*lcoord(2) )   &
137        +( 1.0D0-eta(na)*eta(na) )      &
138        *( 1.0D0-lcoord(2)*lcoord(2) )
139
140      dn_xi(na)  = ( 0.5D0*xi(na)  )            &
141        *( 1.0D0+xi(na) *lcoord(1) ) &
142        +( 0.5D0*xi(na) *lcoord(1) )  &
143        *xi(na)                      &
144        +( 1.0D0-xi(na) *xi(na)  )    &
145        *( -2.0D0*lcoord(1) )
146      dn_eta(na) = ( 0.5D0*eta(na) )            &
147        *( 1.0D0+eta(na)*lcoord(2) ) &
148        +( 0.5D0*eta(na)*lcoord(2) )  &
149        *eta(na)                     &
150        +( 1.0D0-eta(na)*eta(na) )    &
151        *( -2.0D0*lcoord(2) )
152
153    end do
154
155    !--------------------------------------------------------------------
156
157    do na = 1, 9
158
159      func(na, 1) = dn_xi(na)*n_eta(na)
160      func(na, 2) = n_xi(na) *dn_eta(na)
161
162    end do
163
164    !--------------------------------------------------------------------
165
166    return
167
168    !####################################################################
169  end subroutine ShapeDeriv_quad9n
170  !####################################################################
171
172  !####################################################################
173  subroutine NodalNaturalCoord_quad9n(nncoord)
174    !####################################################################
175
176    real(kind = kreal), intent(out) :: nncoord(9, 2)
177
178    !--------------------------------------------------------------------
179
180    ! xi-coordinate at a node in a local element
181    nncoord(1, 1)  = -1.0D0
182    nncoord(2, 1)  =  1.0D0
183    nncoord(3, 1)  =  1.0D0
184    nncoord(4, 1)  = -1.0D0
185    nncoord(5, 1)  =  0.0D0
186    nncoord(6, 1)  =  1.0D0
187    nncoord(7, 1)  =  0.0D0
188    nncoord(8, 1)  = -1.0D0
189    nncoord(9, 1)  =  0.0D0
190    ! eta-coordinate at a node in a local element
191    nncoord(1, 2) = -1.0D0
192    nncoord(2, 2) = -1.0D0
193    nncoord(3, 2) =  1.0D0
194    nncoord(4, 2) =  1.0D0
195    nncoord(5, 2) = -1.0D0
196    nncoord(6, 2) =  0.0D0
197    nncoord(7, 2) =  1.0D0
198    nncoord(8, 2) =  0.0D0
199    nncoord(9, 2) =  0.0D0
200
201    !--------------------------------------------------------------------
202
203    return
204
205    !####################################################################
206  end subroutine NodalNaturalCoord_quad9n
207  !####################################################################
208
209end module shape_quad9n
210