SFEMaNS  version 5.3
Reference documentation for SFEMaNS
GP_2D_p3.f90
Go to the documentation of this file.
1 !===
2 !Author: Jean-Luc Guermond, Copyright December 24th 2018
3 !===
4 MODULE gp_2d_p3
5  PRIVATE
6  PUBLIC element_2d_p3, element_2d_p3_boundary, element_1d_p3
7 CONTAINS
8  SUBROUTINE element_2d_p3 (w, d, p, n_w, l_G)
9  !===triangular element with cubic interpolation
10  !===and seven Gauss integration points
11  !===w(n_w, l_G) : values of shape functions at Gauss points
12  !===d(2, n_w, l_G) : derivatives values of shape functions at Gauss points
13  !===p(l_G) : weight for Gaussian quadrature at Gauss points
14  !=== 3
15  !=== 7 5
16  !=== 6 10 4
17  !=== 1 8 9 2
18  IMPLICIT NONE
19  INTEGER, INTENT(IN) :: n_w, l_G
20  REAL(KIND=8), DIMENSION( n_w, l_G), INTENT(OUT) :: w
21  REAL(KIND=8), DIMENSION(2, n_w, l_G), INTENT(OUT) :: d
22  REAL(KIND=8), DIMENSION(l_G), INTENT(OUT) :: p
23  REAL(KIND=8), DIMENSION(l_G) :: xx, yy
24  INTEGER :: j
25  REAL(KIND=8) :: one=1.d0, two=2.d0, three=3.d0, five=5.d0, nine=9.d0
26  REAL(KIND=8) :: f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, &
27  df1x, df2x, df3x, df4x, df5x, df6x, df7x, df8x, df9x, df10x, &
28  df1y, df2y, df3y, df4y, df5y, df6y, df7y, df8y, df9y, df10y, &
29  x, y, r, a, s1, s2, t1, t2, b1, b2, area, sq
30 
31  f1(x,y) = -0.9d1/0.2d1*x**3 - 0.27d2/0.2d1*x**2*y - 0.27d2/0.2d1*x*y**2 &
32  - 0.9d1/0.2d1*y**3 + 0.9d1*x**2 + 0.18d2*x*y + 0.9d1*y**2 &
33  - 0.11d2/0.2d1*x - 0.11d2/0.2d1*y + 0.1d1
34  f2(x,y) = 0.9d1/0.2d1*x**3 - 0.9d1/0.2d1*x**2 + x
35  f3(x,y) = 0.9d1/0.2d1*y**3 - 0.9d1/0.2d1*y**2 + y
36  f4(x,y) = 0.27d2/0.2d1*x**2*y - 0.9d1/0.2d1*x*y
37  f5(x,y) = 0.27d2/0.2d1*x*y**2 - 0.9d1/0.2d1*x*y
38  f6(x,y) = 0.27d2/0.2d1*x**2*y + 0.27d2*x*y**2 + 0.27d2 / 0.2d1*y**3 &
39  - 0.45d2/0.2d1*x*y - 0.45d2/0.2d1*y**2 + 0.9d1*y
40  f7(x,y) = -0.27d2/0.2d1*x*y**2 - 0.27d2/0.2d1*y**3 + 0.9d1/0.2d1*x*y &
41  + 0.18d2*y**2 - 0.9d1/0.2d1*y
42  f8(x,y) = 0.27d2/0.2d1*x**3 + 0.27d2*x**2*y + 0.27d2/0.2d1*x*y**2 &
43  - 0.45d2/0.2d1*x**2 - 0.45d2/0.2d1*x*y + 0.9d1*x
44  f9(x,y) = -0.27d2/0.2d1*x**3 - 0.27d2/0.2d1*x**2*y + 0.18d2*x**2 &
45  + 0.9d1/0.2d1*x*y - 0.9d1/0.2d1*x
46  f10(x,y) = -27*x**2*y - 27*x*y**2 + 27*x*y
47 
48  df1x(x,y) = -0.27d2/0.2d1*x**2 - 0.27d2*x*y - 0.27d2/0.2d1*y**2 &
49  + 0.18d2*x + 0.18d2*y - 0.11d2/0.2d1
50  df2x(x,y) = 0.27d2/0.2d1*x**2 - 0.9d1*x + 0.1d1
51  df3x(x,y) = 0
52  df4x(x,y) = 27*x*y - 0.9d1/0.2d1 *y
53  df5x(x,y) = 0.27d2/0.2d1*y**2 - 0.9d1/0.2d1*y
54  df6x(x,y) = 27*x*y + 27*y**2 - 0.45d2/0.2d1*y
55  df7x(x,y) = -0.27d2/0.2d1*y**2 + 0.9d1/0.2d1*y
56  df8x(x,y) = 0.81d2/0.2d1*x**2 + 0.54d2*x*y - 0.45d2*x &
57  + 0.27d2/0.2d1*y**2 - 0.45d2/0.2d1*y + 0.9d1
58  df9x(x,y) = -0.81d2/0.2d1*x**2 + 0.36d2*x - 0.27d2*x*y &
59  + 0.9d1/0.2d1*y - 0.9d1/0.2d1
60  df10x(x,y) = -54*x*y - 27*y**2 + 27*y
61 
62  df1y(x,y) = -0.27d2/0.2d1*x**2 - 0.27d2*x*y - 0.27d2/0.2d1* y**2 &
63  + 0.18d2*x + 0.18d2*y - 0.11d2/0.2d1
64  df2y(x,y) = 0.d0
65  df3y(x,y) = 0.27d2/0.2d1*y**2 - 0.9d1*y + 0.1d1
66  df4y(x,y) = 0.27d2/0.2d1*x**2 - 0.9d1/0.2d1*x
67  df5y(x,y) = 27*x*y - 0.9d1/0.2d1*x
68  df6y(x,y) = 54*x*y + 0.81d2/0.2d1*y**2 &
69  - 45*y + 0.27d2/0.2d1*x**2 - 0.45d2/0.2d1*x + 0.9d1
70  df7y(x,y) = -0.81d2/0.2d1*y**2 + 0.36d2*y - 0.27d2*x*y + 0.9d1/0.2d1*x - 0.9d1/0.2d1
71  df8y(x,y) = 27*x**2 + 27*x*y - 0.45d2/0.2d1*x
72  df9y(x,y) = -0.27d2/0.2d1*x**2 + 0.9d1/0.2d1*x
73  df10y(x,y) = -27*x**2 - 54*x*y + 27*x
74 
75  !===Degree 5; 7 Points; Stroud: p. 314, Approximate calculation of
76  !===Multiple integrals (Prentice--Hall), 1971.
77  area = one/two
78  sq = sqrt(three*five)
79  r = one/three; a = area * nine/40
80  s1 = (6 - sq)/21; t1 = (9 + 2*sq)/21; b1 = area * (155 - sq)/1200
81  s2 = (6 + sq)/21; t2 = (9 - 2*sq)/21; b2 = area * (155 + sq)/1200
82  xx(1) = r; yy(1) = r; p(1) = a
83  xx(2) = s1; yy(2) = s1; p(2) = b1
84  xx(3) = s1; yy(3) = t1; p(3) = b1
85  xx(4) = t1; yy(4) = s1; p(4) = b1
86  xx(5) = s2; yy(5) = s2; p(5) = b2
87  xx(6) = s2; yy(6) = t2; p(6) = b2
88  xx(7) = t2; yy(7) = s2; p(7) = b2
89 
90  DO j = 1, l_g
91  w(1, j) = f1(xx(j), yy(j))
92  d(1, 1, j) = df1x(xx(j), yy(j))
93  d(2, 1, j) = df1y(xx(j), yy(j))
94  w(2, j) = f2(xx(j), yy(j))
95  d(1, 2, j) = df2x(xx(j), yy(j))
96  d(2, 2, j) = df2y(xx(j), yy(j))
97  w(3, j) = f3(xx(j), yy(j))
98  d(1, 3, j) = df3x(xx(j), yy(j))
99  d(2, 3, j) = df3y(xx(j), yy(j))
100  w(4, j) = f4(xx(j), yy(j))
101  d(1, 4, j) = df4x(xx(j), yy(j))
102  d(2, 4, j) = df4y(xx(j), yy(j))
103  w(5, j) = f5(xx(j), yy(j))
104  d(1, 5, j) = df5x(xx(j), yy(j))
105  d(2, 5, j) = df5y(xx(j), yy(j))
106  w(6, j) = f6(xx(j), yy(j))
107  d(1, 6, j) = df6x(xx(j), yy(j))
108  d(2, 6, j) = df6y(xx(j), yy(j))
109  w(7, j) = f7(xx(j), yy(j))
110  d(1, 7, j) = df7x(xx(j), yy(j))
111  d(2, 7, j) = df7y(xx(j), yy(j))
112  w(8, j) = f8(xx(j), yy(j))
113  d(1, 8, j) = df8x(xx(j), yy(j))
114  d(2, 8, j) = df8y(xx(j), yy(j))
115  w(9, j) = f9(xx(j), yy(j))
116  d(1, 9, j) = df9x(xx(j), yy(j))
117  d(2, 9, j) = df9y(xx(j), yy(j))
118  w(10, j) = f10(xx(j), yy(j))
119  d(1, 10, j)=df10x(xx(j), yy(j))
120  d(2, 10, j)=df10y(xx(j), yy(j))
121  ENDDO
122 
123  END SUBROUTINE element_2d_p3
124 
125  SUBROUTINE element_2d_p3_boundary (face, d, w, n_w, l_Gs)
126  !===triangular element with cubic interpolation
127  !===and seven Gauss integration points
128  !===d(2, n_w, l_G) : derivatives values of shape functions at Gauss points
129  IMPLICIT NONE
130  INTEGER, INTENT(IN) :: n_w, l_Gs
131  INTEGER, INTENT(IN) :: face
132  REAL(KIND=8), DIMENSION(2, n_w, l_Gs), INTENT(OUT) :: d
133  REAL(KIND=8), DIMENSION( n_w, l_Gs), INTENT(OUT) :: w
134  REAL(KIND=8), DIMENSION(l_Gs) :: xx, yy, gp
135  INTEGER :: j, l
136  REAL(KIND=8) :: half=0.5d0
137  REAL(KIND=8) :: f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, &
138  df1x, df2x, df3x, df4x, df5x, df6x, df7x, df8x, df9x, df10x, &
139  df1y, df2y, df3y, df4y, df5y, df6y, df7y, df8y, df9y, df10y, &
140  x, y
141 
142  f1(x,y) = -0.9d1/0.2d1*x**3 - 0.27d2/0.2d1*x**2*y - 0.27d2/0.2d1*x*y**2 &
143  - 0.9d1/0.2d1*y**3 + 0.9d1*x**2 + 0.18d2*x*y + 0.9d1*y**2 &
144  - 0.11d2/0.2d1*x - 0.11d2/0.2d1*y + 0.1d1
145  f2(x,y) = 0.9d1/0.2d1*x**3 - 0.9d1/0.2d1*x**2 + x
146  f3(x,y) = 0.9d1/0.2d1*y**3 - 0.9d1/0.2d1*y**2 + y
147  f4(x,y) = 0.27d2/0.2d1*x**2*y - 0.9d1/0.2d1*x*y
148  f5(x,y) = 0.27d2/0.2d1*x*y**2 - 0.9d1/0.2d1*x*y
149  f6(x,y) = 0.27d2/0.2d1*x**2*y + 0.27d2*x*y**2 + 0.27d2/0.2d1*y**3 &
150  - 0.45d2/0.2d1*x*y - 0.45d2/0.2d1*y**2 + 0.9d1*y
151  f7(x,y) = -0.27d2/0.2d1*x*y**2 - 0.27d2/0.2d1*y**3 + 0.9d1/0.2d1*x*y &
152  + 0.18d2*y**2 - 0.9d1/0.2d1*y
153  f8(x,y) = 0.27d2/0.2d1*x**3 + 0.27d2*x**2*y + 0.27d2/0.2d1*x*y**2 &
154  - 0.45d2/0.2d1*x**2 - 0.45d2/0.2d1*x*y + 0.9d1*x
155  f9(x,y) = -0.27d2/0.2d1*x**3 - 0.27d2/0.2d1*x**2*y + 0.18d2*x**2 &
156  + 0.9d1/0.2d1*x*y - 0.9d1/0.2d1*x
157  f10(x,y) = -27*x**2*y - 27*x*y**2 + 27*x*y
158 
159  df1x(x,y) = -0.27d2/0.2d1*x**2 - 0.27d2*x*y - 0.27d2/0.2d1*y**2 &
160  + 0.18d2*x + 0.18d2*y - 0.11d2/0.2d1
161  df2x(x,y) = 0.27d2/0.2d1*x**2 - 0.9d1*x + 0.1d1
162  df3x(x,y) = 0
163  df4x(x,y) = 27*x*y - 0.9d1/0.2d1 *y
164  df5x(x,y) = 0.27d2/0.2d1*y**2 - 0.9d1/0.2d1*y
165  df6x(x,y) = 27*x*y + 27*y**2 - 0.45d2/0.2d1*y
166  df7x(x,y) = -0.27d2/0.2d1*y**2 + 0.9d1/0.2d1*y
167  df8x(x,y) = 0.81d2/0.2d1*x**2 + 0.54d2*x*y - 0.45d2*x &
168  + 0.27d2/0.2d1*y**2 - 0.45d2/0.2d1*y + 0.9d1
169  df9x(x,y) = -0.81d2/0.2d1*x**2 + 0.36d2*x - 0.27d2*x*y &
170  + 0.9d1/0.2d1*y - 0.9d1/0.2d1
171  df10x(x,y) = -54*x*y - 27*y**2 + 27*y
172 
173  df1y(x,y) = -0.27d2/0.2d1*x**2 - 0.27d2*x*y - 0.27d2/0.2d1*y**2 &
174  + 0.18d2*x + 0.18d2*y - 0.11d2/0.2d1
175  df2y(x,y) = 0.d0
176  df3y(x,y) = 0.27d2/0.2d1*y**2 - 0.9d1*y + 0.1d1
177  df4y(x,y) = 0.27d2/0.2d1*x**2 - 0.9d1/0.2d1*x
178  df5y(x,y) = 27*x*y - 0.9d1/0.2d1*x
179  df6y(x,y) = 54*x*y + 0.81d2/0.2d1*y**2 &
180  - 45*y + 0.27d2/0.2d1*x**2 - 0.45d2/0.2d1*x + 0.9d1
181  df7y(x,y) = -0.81d2/0.2d1*y**2 + 0.36d2*y - 0.27d2*x*y + 0.9d1/0.2d1*x - 0.9d1/0.2d1
182  df8y(x,y) = 27*x**2 + 27*x*y - 0.45d2/0.2d1*x
183  df9y(x,y) = -0.27d2/0.2d1*x**2 + 0.9d1/0.2d1*x
184  df10y(x,y) = -27*x**2 - 54*x*y + 27*x
185 
186  gp(1) = -sqrt((15.d0+2*sqrt(30.d0))/35.d0)
187  gp(2) = -sqrt((15.d0-2*sqrt(30.d0))/35.d0)
188  gp(3) = sqrt((15.d0-2*sqrt(30.d0))/35.d0)
189  gp(4) = sqrt((15.d0+2*sqrt(30.d0))/35.d0)
190 
191  IF (face==1) THEN
192  DO l = 1, l_gs
193  xx(l) = 1.d0-(gp(l)+1.d0)*half
194  yy(l) = 0.d0+(gp(l)+1.d0)*half
195  END DO
196  ELSE IF (face==2) THEN
197  DO l = 1, l_gs
198  xx(l) = 0.d0
199  yy(l) = (gp(l)+1.d0)*half
200  END DO
201  ELSE
202  DO l = 1, l_gs
203  xx(l) = (gp(l)+1.d0)*half
204  yy(l) = 0.d0
205  END DO
206  END IF
207 
208  DO j = 1, l_gs
209  w(1, j) = f1(xx(j), yy(j))
210  d(1, 1, j) = df1x(xx(j), yy(j))
211  d(2, 1, j) = df1y(xx(j), yy(j))
212  w(2, j) = f2(xx(j), yy(j))
213  d(1, 2, j) = df2x(xx(j), yy(j))
214  d(2, 2, j) = df2y(xx(j), yy(j))
215  w(3, j) = f3(xx(j), yy(j))
216  d(1, 3, j) = df3x(xx(j), yy(j))
217  d(2, 3, j) = df3y(xx(j), yy(j))
218  w(4, j) = f4(xx(j), yy(j))
219  d(1, 4, j) = df4x(xx(j), yy(j))
220  d(2, 4, j) = df4y(xx(j), yy(j))
221  w(5, j) = f5(xx(j), yy(j))
222  d(1, 5, j) = df5x(xx(j), yy(j))
223  d(2, 5, j) = df5y(xx(j), yy(j))
224  w(6, j) = f6(xx(j), yy(j))
225  d(1, 6, j) = df6x(xx(j), yy(j))
226  d(2, 6, j) = df6y(xx(j), yy(j))
227  w(7, j) = f7(xx(j), yy(j))
228  d(1, 7, j) = df7x(xx(j), yy(j))
229  d(2, 7, j) = df7y(xx(j), yy(j))
230  w(8, j) = f8(xx(j), yy(j))
231  d(1, 8, j) = df8x(xx(j), yy(j))
232  d(2, 8, j) = df8y(xx(j), yy(j))
233  w(9, j) = f9(xx(j), yy(j))
234  d(1, 9, j) = df9x(xx(j), yy(j))
235  d(2, 9, j) = df9y(xx(j), yy(j))
236  w(10, j) = f10(xx(j), yy(j))
237  d(1, 10, j)=df10x(xx(j), yy(j))
238  d(2, 10, j)=df10y(xx(j), yy(j))
239  ENDDO
240 
241  END SUBROUTINE element_2d_p3_boundary
242 
243  SUBROUTINE element_1d_p3(w, d, p, n_ws, l_Gs)
244  !===one-dimensional element with cubic interpolation
245  !===and 4 Gauss integration points
246  !===w(n_w, l_G) : values of shape functions at Gauss points
247  !===d(1, n_w, l_G) : derivatives values of shape functions at Gauss points
248  !===p(l_G) : weight for Gaussian quadrature at Gauss points
249  !===Enumeration: 1 3 4 2
250  IMPLICIT NONE
251  INTEGER, INTENT(IN) :: n_ws, l_Gs
252  REAL(KIND=8), DIMENSION( n_ws, l_Gs), INTENT(OUT) :: w
253  REAL(KIND=8), DIMENSION(1, n_ws, l_Gs), INTENT(OUT) :: d
254  REAL(KIND=8), DIMENSION(l_Gs), INTENT(OUT) :: p
255  REAL(KIND=8), DIMENSION(l_Gs) :: xx
256  INTEGER :: j
257  REAL(KIND=8) :: f1, f2, f3, f4, df1, df2, df3, df4, x
258 
259  f1(x) = -0.1d1/0.16d2 + x/0.16d2 + 0.9d1/0.16d2*x**2 - 0.9d1/0.16d2*x**3
260  f2(x) = -x/0.16d2 + 0.9d1/0.16d2*x**2 + 0.9d1/0.16d2*x**3 - 0.1d1/0.16d2
261  f3(x) = -0.9d1/0.16d2*x**2 + 0.27d2/0.16d2*x**3 + 0.9d1/0.16d2 - 0.27d2/0.16d2*x
262  f4(x) = -0.9d1/0.16d2*x**2 - 0.27d2/0.16d2*x**3 + 0.9d1/0.16d2 + 0.27d2/0.16d2*x
263 
264  df1(x) = 0.1d1/0.16d2 - 0.27d2/0.16d2*x**2 + 0.9d1/0.8d1*x
265  df2(x) = -0.1d1/0.16d2 + 0.27d2/0.16d2*x**2 + 0.9d1/0.8d1*x
266  df3(x) = -0.27d2/0.16d2 - 0.9d1/0.8d1*x + 0.81d2/0.16d2*x**2
267  df4(x) = -0.9d1/0.8d1*x - 0.81d2/0.16d2*x**2 + 0.27d2/0.16d2
268 
269  xx(1) = -sqrt((15.d0+2*sqrt(30.d0))/35.d0)
270  xx(2) = -sqrt((15.d0-2*sqrt(30.d0))/35.d0)
271  xx(3) = sqrt((15.d0-2*sqrt(30.d0))/35.d0)
272  xx(4) = sqrt((15.d0+2*sqrt(30.d0))/35.d0)
273  p(1) = (18.d0-sqrt(30.d0))/36.d0
274  p(2) = (18.d0+sqrt(30.d0))/36.d0
275  p(3) = (18.d0+sqrt(30.d0))/36.d0
276  p(4) = (18.d0-sqrt(30.d0))/36.d0
277 
278  DO j = 1, l_gs
279  w(1, j) = f1(xx(j))
280  d(1, 1, j) = df1(xx(j))
281  w(2, j) = f2(xx(j))
282  d(1, 2, j) = df2(xx(j))
283  w(3, j) = f3(xx(j))
284  d(1, 3, j) = df3(xx(j))
285  w(4, j) = f4(xx(j))
286  d(1, 4, j) = df4(xx(j))
287  ENDDO
288  END SUBROUTINE element_1d_p3
289 END MODULE gp_2d_p3
290 
gp_2d_p3::element_2d_p3_boundary
subroutine, public element_2d_p3_boundary(face, d, w, n_w, l_Gs)
Definition: GP_2D_p3.f90:126
gp_2d_p3
Definition: GP_2D_p3.f90:4
gp_2d_p3::element_1d_p3
subroutine, public element_1d_p3(w, d, p, n_ws, l_Gs)
Definition: GP_2D_p3.f90:244
gp_2d_p3::element_2d_p3
subroutine, public element_2d_p3(w, d, p, n_w, l_G)
Definition: GP_2D_p3.f90:9