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
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
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
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
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
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
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
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
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))
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
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, &
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
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
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
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
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
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)
193 xx(l) = 1.d0-(gp(l)+1.d0)*half
194 yy(l) = 0.d0+(gp(l)+1.d0)*half
196 ELSE IF (face==2)
THEN 199 yy(l) = (gp(l)+1.d0)*half
203 xx(l) = (gp(l)+1.d0)*half
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))
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
257 REAL(KIND=8) :: f1, f2, f3, f4, df1, df2, df3, df4, x
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
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
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
280 d(1, 1, j) = df1(xx(j))
282 d(1, 2, j) = df2(xx(j))
284 d(1, 3, j) = df3(xx(j))
286 d(1, 4, j) = df4(xx(j))
subroutine, public element_2d_p3_boundary(face, d, w, n_w, l_Gs)
subroutine, public element_1d_p3(w, d, p, n_ws, l_Gs)
subroutine, public element_2d_p3(w, d, p, n_w, l_G)