SFEMaNS  version 4.1 (work in progress)
Reference documentation for SFEMaNS
 All Classes Files Functions Variables Groups Pages
gauss_points_2d_p2.f90
Go to the documentation of this file.
1 !
2 !Authors: Jean-Luc Guermond, Lugi Quartapelle, Copyright 1994
3 !
4 MODULE mod_gauss_points_2d_p2
5 
6  PRIVATE
7 
8  PUBLIC gauss_points_2d_p2
9 
10 CONTAINS
11 
12  SUBROUTINE gauss_points_2d_p2(mesh)
13 
14  USE def_type_mesh
15 
16  IMPLICIT NONE
17 
18  TYPE(mesh_type), TARGET :: mesh
19 
20  INTEGER, POINTER :: me, mes
21  INTEGER, DIMENSION(:, :), POINTER :: jj
22  INTEGER, DIMENSION(:, :), POINTER :: js
23  REAL(KIND=8), DIMENSION(:, :), POINTER :: rr
24 
25  INTEGER, PARAMETER :: k_d = 2, n_w = 6, l_g = 7, &
26  n_ws = 3, l_gs = 3
27 
28  REAL(KIND=8), DIMENSION(:, :), POINTER :: ww
29  REAL(KIND=8), DIMENSION(:, :), POINTER :: wws
30  REAL(KIND=8), DIMENSION(:, :, :, :), POINTER :: dw
31  REAL(KIND=8), DIMENSION(:, :, :), POINTER :: rnorms
32  REAL(KIND=8), DIMENSION(:, :), POINTER :: rj
33  REAL(KIND=8), DIMENSION(:, :), POINTER :: rjs
34  REAL(KIND=8), DIMENSION(:, :, :, :), POINTER :: dw_s
35  REAL(KIND=8), DIMENSION(:, :, :, :), POINTER :: dwps !special!
36  REAL(KIND=8), DIMENSION(:, :, :, :), POINTER :: dws !SPECIAL!
37 
38  REAL(KIND=8), DIMENSION(k_d, n_w, l_G) :: dd
39  REAL(KIND=8), DIMENSION(1, n_ws, l_Gs) :: dds
40  REAL(KIND=8), DIMENSION(l_G) :: pp
41  REAL(KIND=8), DIMENSION(l_Gs) :: pps
42  REAL(KIND=8), DIMENSION(n_w) :: r
43  REAL(KIND=8), DIMENSION(n_ws) :: rs
44  REAL(KIND=8), DIMENSION(k_d, k_d) :: dr
45  REAL(KIND=8), DIMENSION(1, k_d) :: drs
46 
47  REAL(KIND=8) :: rjac, rjacs, x
48  INTEGER :: m, l, k, k1, n, ms, ls, face, n1, n2, cote, sign
49  REAL(KIND=8), DIMENSION(2) :: rnor, rsd ! JLG April 2009
50 
51  !TEST
52  REAL(KIND=8), DIMENSION(n_w, l_G) :: ww_s
53  !TEST
54  !NEW Mai 26, 2004
55  INTEGER, DIMENSION(n_ws) :: nface
56  INTEGER, DIMENSION(l_Gs) :: ngauss
57  !NEW Mai 26, 2004
58  me => mesh%me
59  mes => mesh%mes
60  jj => mesh%jj
61  js => mesh%jjs
62  rr => mesh%rr
63 
64  ! JLG (April 2009)
65  IF (mesh%edge_stab) THEN
66  ALLOCATE(mesh%gauss%dwni(n_w, l_gs, 2, mesh%mi))
67  ALLOCATE(mesh%gauss%rji(l_gs, mesh%mi))
68  ALLOCATE(mesh%gauss%rnormsi(k_d, l_gs, 2, mesh%mi)) !JLG June 4 2012
69  ALLOCATE(mesh%gauss%wwsi(n_ws, l_gs)) !JLG June 4 2012
70  END IF
71  ! JLG (April 2009)
72 
73  ALLOCATE(mesh%gauss%ww(n_w, l_g))
74  ALLOCATE(mesh%gauss%wws(n_ws, l_gs))
75  ALLOCATE(mesh%gauss%dw(k_d, n_w, l_g, me ))
76  ALLOCATE(mesh%gauss%rnorms(k_d, l_gs, mes))
77  ALLOCATE(mesh%gauss%rj(l_g, me ))
78  ALLOCATE(mesh%gauss%rjs(l_gs, mes))
79  !NEW Juin 8, 2004
80  ALLOCATE(mesh%gauss%dw_s(k_d, n_w, l_gs, mesh%mes))
81  !NEW Juin 8, 2004
82  ALLOCATE(mesh%gauss%dwps(1, n_ws, l_gs, mes))
83  ALLOCATE(mesh%gauss%dws(1, n_ws, l_gs, mes))
84 
85  ww => mesh%gauss%ww
86  wws => mesh%gauss%wws
87  dw => mesh%gauss%dw
88  rnorms => mesh%gauss%rnorms
89  rj => mesh%gauss%rj
90  rjs => mesh%gauss%rjs
91  !NEW Mai 26, 2004
92  dw_s => mesh%gauss%dw_s
93  !NEW Mai 26, 2004
94  dwps => mesh%gauss%dwps
95  dws => mesh%gauss%dws
96 
97  mesh%gauss%k_d = k_d
98  mesh%gauss%n_w = n_w
99  mesh%gauss%l_G = l_g
100  mesh%gauss%n_ws = n_ws
101  mesh%gauss%l_Gs = l_gs
102 
103  ! evaluate and store the values of derivatives and of the
104  ! jacobian determinant at Gauss points of all volume elements
105 
106  ! volume elements
107 
108  CALL element_2d_p2(ww, dd, pp)
109 
110  DO m = 1, me
111 
112  DO l = 1, l_g
113 
114  DO k = 1, k_d
115  r = rr(k, jj(:,m))
116  DO k1 = 1, k_d
117  dr(k, k1) = sum(r * dd(k1,:,l))
118  ENDDO
119  ENDDO
120 
121  rjac = dr(1,1)*dr(2,2) - dr(1,2)*dr(2,1)
122 
123  DO n = 1, n_w
124  dw(1, n, l, m) &
125  = (+ dd(1,n,l)*dr(2,2) - dd(2,n,l)*dr(2,1))/rjac
126  dw(2, n, l, m) &
127  = (- dd(1,n,l)*dr(1,2) + dd(2,n,l)*dr(1,1))/rjac
128  ENDDO
129 
130  rj(l, m) = rjac * pp(l)
131 
132  ENDDO
133 
134  ENDDO
135 
136  !TEST
137  CALL element_1d_p2(wws, dds, pps)
138  !TEST
139 
140  ! surface elements
141  DO ms = 1, mes
142  m = mesh%neighs(ms)
143  !Il faut savoir quelle face est la bonne
144  DO n = 1, 3
145  IF (minval(abs(js(:,ms)-jj(n,m)))==0) cycle
146  face = n
147  END DO
148  CALL element_2d_p2_boundary(face, dd, ww_s)
149  n1 = modulo(face,3)+1; n2 = modulo(face+1,3)+1
150  IF (js(1,ms)==jj(n1,m) .AND. js(2,ms)==jj(n2,m)) THEN
151  nface(1)=n1; nface(2)=n2; nface(3)=face+3
152  ngauss(1)=1; ngauss(2)=2; ngauss(3)=3;
153  sign = 1
154  ELSE IF (js(1,ms)==jj(n2,m) .AND. js(2,ms)==jj(n1,m)) THEN
155  nface(1)=n2; nface(2)=n1; nface(3)=face+3
156  ngauss(1)=3; ngauss(2)=2; ngauss(3)=1;
157  sign = -1
158  ELSE
159  WRITE(*,*) ' BUG : gauss_points_2d_p2'
160  stop
161  END IF
162 
163  DO ls = 1, l_gs
164 
165  !write(*,*) wws(:,ls)
166  !write(*,*) ww_s(nface,ngauss(ls)), sign
167 
168  DO k = 1, k_d
169  r = rr(k, jj(:,m))
170  DO k1 = 1, k_d
171  dr(k, k1) = sum(r * dd(k1,:,ls))
172  ENDDO
173  ENDDO
174 
175  rjac = (dr(1,1)*dr(2,2) - dr(1,2)*dr(2,1))
176 
177  DO n = 1, n_w
178  dw_s(1, n, ngauss(ls), ms) &
179  = (+ dd(1,n,ls)*dr(2,2) - dd(2,n,ls)*dr(2,1))/rjac
180  dw_s(2, n, ngauss(ls), ms) &
181  = (- dd(1,n,ls)*dr(1,2) + dd(2,n,ls)*dr(1,1))/rjac
182  ENDDO
183 
184  ENDDO
185 
186  ENDDO
187 
188  CALL element_1d_p2(wws, dds, pps)
189 
190  DO ms = 1, mes
191 
192  DO ls = 1, l_gs
193 
194 
195  DO k = 1, k_d
196  rs = rr(k, js(:,ms))
197  drs(1, k) = sum(rs * dds(1,:,ls))
198  ENDDO
199 
200  rjacs = sqrt( drs(1,1)**2 + drs(1,2)**2 )
201 
202  rnorms(1, ls, ms) = - drs(1,2)/rjacs ! outward normal
203  rnorms(2, ls, ms) = + drs(1,1)/rjacs ! outward normal
204  !erreur de signe corrigee le 10 Oct 2003. La normale etait interieure
205  rjs(ls, ms) = rjacs * pps(ls)
206 
207  m = mesh%neighs(ms)
208  !Il faut savoir quelle face est la bonne
209  DO n = 1, n_w
210  IF (minval(abs(js(:,ms)-jj(n,m)))==0) cycle
211  face = n
212  END DO
213  rs(1:k_d) = rr(:,jj(face,m)) - (rr(:,js(1,ms))+rr(:,js(2,ms)))/2
214  x = sum(rnorms(:,ls,ms)*rs(1:k_d))
215  IF (x>0) THEN
216  rnorms(:,ls,ms) = -rnorms(:,ls,ms)
217  !write(*,*) 'retournement de n', ms, ls
218  !erreur de signe corrigee le 8 juin 2004
219  END IF
220 
221  ENDDO
222 
223  ENDDO
224 
225 
226  DO ms = 1, mes ! necessary only for evaluating gradient
227  ! tangential to the surface (ex COMMON Gauss_tan)
228 
229  DO ls = 1, l_gs
230  !TEST
231  !dwps(1, :, ls, ms) = dds(1, :, ls) / (rjs(ls, ms)/pps(ls))
232  !do n = 1, n_ws
233  !rjac = SUM(rnorms(:,ls,ms)*dw_s(:,n,ls,ms))
234  !write(*,*) n, ' js ', js(:,ms), rjs(ls, ms)/pps(ls)
235  !write(*,*) ms, ls, dw_s(1, n, ls, ms) - rjac*rnorms(1,ls,ms), dw_s(2, n, ls, ms) - rjac*rnorms(2,ls,ms)
236  !write(*,*) ms, ls, dwps(1, n, ls, ms)*rnorms(2,ls,ms), -dwps(1, n, ls, ms)*rnorms(1,ls,ms), &
237  ! SQRT(SUM(rnorms(:,ls,ms)**2)),rnorms(:,ls,ms)
238  !write(*,*) ' '
239  !end do
240  !TEST
241  dwps(1, :, ls, ms) = dds(1, :, ls) * pps(ls)
242  dws(1, :, ls, ms) = dds(1, :, ls)
243 
244  ENDDO
245 
246  ENDDO
247 
248 
249 
250  !-----------------------------------------------------------------------------
251  ! CELL INTERFACES (JLG, April 2009)
252  IF (mesh%edge_stab) THEN
253  CALL element_1d_p2(mesh%gauss%wwsi, dds, pps) !JLG June 4 2012
254  DO ms = 1, mesh%mi
255  DO cote = 1, 2
256  m = mesh%neighi(cote,ms)
257  !Il faut savoir quelle face est la bonne
258  DO n = 1, 3
259  IF (minval(abs(mesh%jjsi(:,ms)-jj(n,m)))==0) cycle
260  face = n
261  END DO
262  CALL element_2d_p2_boundary(face, dd, ww_s)
263  n1 = modulo(face,3)+1; n2 = modulo(face+1,3)+1
264  IF (mesh%jjsi(1,ms)==jj(n1,m) .AND. mesh%jjsi(2,ms)==jj(n2,m)) THEN
265  nface(1)=n1; nface(2)=n2
266  !ngauss(1)=1; ngauss(2)=2 ! JLG: BUG Jan 26 2010
267  ngauss(1)=1; ngauss(2)=2; ngauss(3)=3
268  ELSE IF (mesh%jjsi(1,ms)==jj(n2,m) .AND. mesh%jjsi(2,ms)==jj(n1,m)) THEN
269  nface(1)=n2; nface(2)=n1
270  !ngauss(1)=2; ngauss(2)=1 ! JLG: BUG Jan 26 2010
271  ngauss(1)=3; ngauss(2)=2; ngauss(3)=1
272  ELSE
273  WRITE(*,*) ' BUG : gauss_points_2d_p2'
274  stop
275  END IF
276  DO ls = 1, l_gs
277  ! Compute normal vector
278  DO k = 1, k_d
279  rs = rr(k, mesh%jjsi(:,ms))
280  drs(1, k) = sum(rs * dds(1,:,ls))
281  ENDDO
282  rjacs = sqrt( drs(1,1)**2 + drs(1,2)**2 )
283  mesh%gauss%rji(ls, ms) = rjacs * pps(ls)
284  rnor(1) = drs(1,2)/rjacs
285  rnor(2) = - drs(1,1)/rjacs
286  mesh%gauss%rnormsi(1, ls, cote, ms) = rnor(1) !JLG, June 4 2012
287  mesh%gauss%rnormsi(2, ls, cote, ms) = rnor(2) !JLG, June 4 2012
288  rsd = rr(:,jj(face,m)) - (rr(:,mesh%jjsi(1,ms))+rr(:,mesh%jjsi(2,ms)))/2
289  x = sum(rnor(:)*rsd)
290  IF (x>0) THEN ! Verify the normal is outward
291  rnor = -rnor
292  !WRITE(*,*) 'retournement de n', ms, ls, cote, x
293  END IF
294  ! End computation normal vector
295  !WRITE(*,*) SUM((rr(:,mesh%jjsi(1,ms))-rr(:,mesh%jjsi(2,ms)))*rnor)
296  DO k = 1, k_d
297  r = rr(k, jj(:,m))
298  DO k1 = 1, k_d
299  dr(k, k1) = sum(r * dd(k1,:,ls))
300  ENDDO
301  ENDDO
302  rjac = (dr(1,1)*dr(2,2) - dr(1,2)*dr(2,1))
303  DO n = 1, n_w
304  !mesh%gauss%dwi(1, n, ngauss(ls), cote, ms) &
305  ! = (+ dd(1,n,ls)*dr(2,2) - dd(2,n,ls)*dr(2,1))/rjac
306  !mesh%gauss%dwi(2, n, ngauss(ls), cote, ms) &
307  ! = (- dd(1,n,ls)*dr(1,2) + dd(2,n,ls)*dr(1,1))/rjac
308  !r(n) =(dd(1,n,ls)*dr(2,2) - dd(2,n,ls)*dr(2,1))/rjac
309 
310  mesh%gauss%dwni(n, ngauss(ls), cote, ms) &
311  = rnor(1)*(+ dd(1,n,ls)*dr(2,2) - dd(2,n,ls)*dr(2,1))/rjac &
312  + rnor(2)*(- dd(1,n,ls)*dr(1,2) + dd(2,n,ls)*dr(1,1))/rjac
313  ENDDO
314  !IF (ABS(r(1) + 0.5*(r(5)+r(6))- (dw(1,1,ls,m)+0.5*(dw(1,5,ls,m)+dw(1,6,ls,m)))).GE.1.d-12 .OR. &
315  ! ABS(r(2) + 0.5*(r(4)+r(6))- (dw(1,2,ls,m)+0.5*(dw(1,4,ls,m)+dw(1,6,ls,m)))).GE.1.d-12 .OR. &
316  ! ABS(r(3) + 0.5*(r(4)+r(5))- (dw(1,3,ls,m)+0.5*(dw(1,4,ls,m)+dw(1,5,ls,m)))).GE.1.d-12) THEN
317  ! write(*,*) ' BUG '
318  !END IF
319  ENDDO
320  END DO
321  ENDDO
322 
323  DO ms = 1, mesh%mi
324  DO cote = 1, 2
325  m = mesh%neighi(cote,ms)
326  DO ls = 1, l_gs
327  x = abs(sqrt(sum(mesh%gauss%rnormsi(:, ls, cote, ms)**2))-1.d0)
328  !write(*,*) ' x ',x, ABS(SQRT(SUM(mesh%gauss%rnormsi(:, ls, cote, ms)**2)))
329  IF (x.GE.1.d-13) THEN
330  write(*,*) 'Bug1 in normal '
331  stop
332  END IF
333  rjac = (rr(1, mesh%jjsi(2,ms)) - rr(1, mesh%jjsi(1,ms)))*mesh%gauss%rnormsi(1, ls, cote, ms) &
334  + (rr(2, mesh%jjsi(2,ms)) - rr(2, mesh%jjsi(1,ms)))*mesh%gauss%rnormsi(2, ls, cote, ms)
335  x = sqrt((rr(1, mesh%jjsi(2,ms)) - rr(1, mesh%jjsi(1,ms)))**2 +&
336  (rr(2, mesh%jjsi(2,ms)) - rr(2, mesh%jjsi(1,ms)))**2)
337  !write(*,*) ' rjac/x ', rjac/x, x, rjac
338  IF (abs(rjac/x) .GE. 1.d-13) THEN
339  write(*,*) 'Bug2 in normal '
340  stop
341  END IF
342  END DO
343  END DO
344  END DO
345  END IF
346  !-----------------------------------------------------------------------------
347  !PRINT*, 'end of gen_Gauss'
348 
349  CONTAINS
350 
351  SUBROUTINE element_2d_p2 (w, d, p)
352 
353  ! triangular element with quadratic interpolation
354  ! and seven Gauss integration points
355 
356  ! w(n_w, l_G) : values of shape functions at Gauss points
357  ! d(2, n_w, l_G) : derivatives values of shape functions at Gauss points
358  ! p(l_G) : weight for Gaussian quadrature at Gauss points
359 
360  IMPLICIT NONE
361 
362  REAL(KIND=8), DIMENSION( n_w, l_G), INTENT(OUT) :: w
363  REAL(KIND=8), DIMENSION(2, n_w, l_G), INTENT(OUT) :: d
364  REAL(KIND=8), DIMENSION(l_G), INTENT(OUT) :: p
365 
366  REAL(KIND=8), DIMENSION(l_G) :: xx, yy
367  INTEGER :: j
368 
369  REAL(KIND=8) :: zero = 0, half = 0.5, one = 1, &
370  two = 2, three = 3, four = 4, &
371  five = 5, nine = 9
372 
373  REAL(KIND=8) :: f1, f2, f3, f4, f5, f6, &
374  df1x, df2x, df3x, df4x, df5x, df6x, &
375  df1y, df2y, df3y, df4y, df5y, df6y, &
376  x, y, r, a, s1, s2, t1, t2, b1, b2, area, sq
377 
378 
379  f1(x, y) = (half - x - y) * (one - x - y) * two
380  f2(x, y) = x * (x - half) * two
381  f3(x, y) = y * (y - half) * two
382  f4(x, y) = x * y * four
383  f5(x, y) = y * (one - x - y) * four
384  f6(x, y) = x * (one - x - y) * four
385 
386  df1x(x, y) = -three + four * (x + y)
387  df2x(x, y) = (two*x - half) * two
388  df3x(x, y) = zero
389  df4x(x, y) = y * four
390  df5x(x, y) = -y * four
391  df6x(x, y) = (one - two*x - y) * four
392 
393  df1y(x, y) = -three + four * (x + y)
394  df2y(x, y) = zero
395  df3y(x, y) = (two*y - half) * two
396  df4y(x, y) = x * four
397  df5y(x, y) = (one - x - two*y) * four
398  df6y(x, y) = -x * four
399 
400  ! Degree 5; 7 Points; Stroud: p. 314, Approximate calculation of
401  ! Multiple integrals (Prentice--Hall), 1971.
402 
403  area = one/two
404 
405  sq = sqrt(three*five)
406 
407  r = one/three; a = area * nine/40
408 
409  s1 = (6 - sq)/21; t1 = (9 + 2*sq)/21; b1 = area * (155 - sq)/1200
410 
411  s2 = (6 + sq)/21; t2 = (9 - 2*sq)/21; b2 = area * (155 + sq)/1200
412 
413  xx(1) = r; yy(1) = r; p(1) = a
414 
415  xx(2) = s1; yy(2) = s1; p(2) = b1
416  xx(3) = s1; yy(3) = t1; p(3) = b1
417  xx(4) = t1; yy(4) = s1; p(4) = b1
418 
419  xx(5) = s2; yy(5) = s2; p(5) = b2
420  xx(6) = s2; yy(6) = t2; p(6) = b2
421  xx(7) = t2; yy(7) = s2; p(7) = b2
422 
423 
424  DO j = 1, l_g
425 
426  w(1, j) = f1(xx(j), yy(j))
427  d(1, 1, j) = df1x(xx(j), yy(j))
428  d(2, 1, j) = df1y(xx(j), yy(j))
429 
430  w(2, j) = f2(xx(j), yy(j))
431  d(1, 2, j) = df2x(xx(j), yy(j))
432  d(2, 2, j) = df2y(xx(j), yy(j))
433 
434  w(3, j) = f3(xx(j), yy(j))
435  d(1, 3, j) = df3x(xx(j), yy(j))
436  d(2, 3, j) = df3y(xx(j), yy(j))
437 
438  w(4, j) = f4(xx(j), yy(j))
439  d(1, 4, j) = df4x(xx(j), yy(j))
440  d(2, 4, j) = df4y(xx(j), yy(j))
441 
442  w(5, j) = f5(xx(j), yy(j))
443  d(1, 5, j) = df5x(xx(j), yy(j))
444  d(2, 5, j) = df5y(xx(j), yy(j))
445 
446  w(6, j) = f6(xx(j), yy(j))
447  d(1, 6, j) = df6x(xx(j), yy(j))
448  d(2, 6, j) = df6y(xx(j), yy(j))
449 
450  ENDDO
451 
452  END SUBROUTINE element_2d_p2
453 
454  SUBROUTINE element_2d_p2_boundary (face, d, w)
455 
456  ! triangular element with quadratic interpolation
457  ! and seven Gauss integration points
458 
459  ! d(2, n_w, l_G) : derivatives values of shape functions at Gauss points
460 
461  IMPLICIT NONE
462 
463  INTEGER, INTENT(IN) :: face
464  REAL(KIND=8), DIMENSION(2, n_w, l_Gs), INTENT(OUT) :: d
465  REAL(KIND=8), DIMENSION( n_w, l_G), INTENT(OUT) :: w
466 
467  REAL(KIND=8), DIMENSION(l_G) :: xx, yy
468  INTEGER :: j
469 
470  REAL(KIND=8) :: zero = 0, half = 0.5d0, one = 1, &
471  two = 2, three = 3, four = 4, &
472  five = 5
473 
474  REAL(KIND=8) :: f1, f2, f3, f4, f5, f6, &
475  df1x, df2x, df3x, df4x, df5x, df6x, &
476  df1y, df2y, df3y, df4y, df5y, df6y, &
477  x, y
478 
479 
480  f1(x, y) = (half - x - y) * (one - x - y) * two
481  f2(x, y) = x * (x - half) * two
482  f3(x, y) = y * (y - half) * two
483  f4(x, y) = x * y * four
484  f5(x, y) = y * (one - x - y) * four
485  f6(x, y) = x * (one - x - y) * four
486 
487  df1x(x, y) = -three + four * (x + y)
488  df2x(x, y) = (two*x - half) * two
489  df3x(x, y) = zero
490  df4x(x, y) = y * four
491  df5x(x, y) = -y * four
492  df6x(x, y) = (one - two*x - y) * four
493 
494  df1y(x, y) = -three + four * (x + y)
495  df2y(x, y) = zero
496  df3y(x, y) = (two*y - half) * two
497  df4y(x, y) = x * four
498  df5y(x, y) = (one - x - two*y) * four
499  df6y(x, y) = -x * four
500 
501  IF (face==1) THEN
502  xx(1) = half + half*sqrt(three/five)
503  yy(1) = half - half*sqrt(three/five)
504 
505  xx(2) = half
506  yy(2) = half
507 
508  xx(3) = half - half*sqrt(three/five)
509  yy(3) = half + half*sqrt(three/five)
510  ELSE IF (face==2) THEN
511  xx(1) = zero
512  yy(1) = half + half*sqrt(three/five)
513 
514  xx(2) = zero
515  yy(2) = half
516 
517  xx(3) = zero
518  yy(3) = half - half*sqrt(three/five)
519  ELSE
520  xx(1) = half - half*sqrt(three/five)
521  yy(1) = zero
522 
523  xx(2) = half
524  yy(2) = zero
525 
526  xx(3) = half + half*sqrt(three/five)
527  yy(3) = zero
528  END IF
529 
530 
531  DO j = 1, l_gs
532 
533  w(1, j) = f1(xx(j), yy(j))
534  d(1, 1, j) = df1x(xx(j), yy(j))
535  d(2, 1, j) = df1y(xx(j), yy(j))
536 
537  w(2, j) = f2(xx(j), yy(j))
538  d(1, 2, j) = df2x(xx(j), yy(j))
539  d(2, 2, j) = df2y(xx(j), yy(j))
540 
541  w(3, j) = f3(xx(j), yy(j))
542  d(1, 3, j) = df3x(xx(j), yy(j))
543  d(2, 3, j) = df3y(xx(j), yy(j))
544 
545  w(4, j) = f4(xx(j), yy(j))
546  d(1, 4, j) = df4x(xx(j), yy(j))
547  d(2, 4, j) = df4y(xx(j), yy(j))
548 
549  w(5, j) = f5(xx(j), yy(j))
550  d(1, 5, j) = df5x(xx(j), yy(j))
551  d(2, 5, j) = df5y(xx(j), yy(j))
552 
553  w(6, j) = f6(xx(j), yy(j))
554  d(1, 6, j) = df6x(xx(j), yy(j))
555  d(2, 6, j) = df6y(xx(j), yy(j))
556 
557  ENDDO
558 
559  END SUBROUTINE element_2d_p2_boundary
560 
561  !------------------------------------------------------------------------------
562 
563  SUBROUTINE element_1d_p2(w, d, p)
564 
565  ! one-dimensional element with linear interpolation
566  ! and two Gauss integration points
567 
568  ! w(n_w, l_G) : values of shape functions at Gauss points
569  ! d(1, n_w, l_G) : derivatives values of shape functions at Gauss points
570  ! p(l_G) : weight for Gaussian quadrature at Gauss points
571 
572  IMPLICIT NONE
573 
574  REAL(KIND=8), DIMENSION( n_ws, l_Gs), INTENT(OUT) :: w
575  REAL(KIND=8), DIMENSION(1, n_ws, l_Gs), INTENT(OUT) :: d
576  REAL(KIND=8), DIMENSION(l_Gs), INTENT(OUT) :: p
577 
578  REAL(KIND=8), DIMENSION(l_Gs) :: xx
579  INTEGER :: j
580 
581  REAL(KIND=8) :: zero = 0, one = 1, two = 2, three = 3, &
582  five = 5, eight= 8, nine = 9
583 
584  REAL(KIND=8) :: f1, f2, f3, df1, df2, df3, x
585 
586  f1(x) = (x - one)*x/two
587  f2(x) = (x + one)*x/two
588  f3(x) = (x + one)*(one - x)
589 
590  df1(x) = (two*x - one)/two
591  df2(x) = (two*x + one)/two
592  df3(x) = -two*x
593 
594  xx(1) = -sqrt(three/five)
595  xx(2) = zero
596  xx(3) = sqrt(three/five)
597 
598  p(1) = five/nine
599  p(2) = eight/nine
600  p(3) = five/nine
601 
602  DO j = 1, l_gs
603 
604  w(1, j) = f1(xx(j))
605  d(1, 1, j) = df1(xx(j))
606 
607  w(2, j) = f2(xx(j))
608  d(1, 2, j) = df2(xx(j))
609 
610  w(3, j) = f3(xx(j))
611  d(1, 3, j) = df3(xx(j))
612 
613 
614  ENDDO
615 
616  END SUBROUTINE element_1d_p2
617 
618  END SUBROUTINE gauss_points_2d_p2
619 
620 END MODULE mod_gauss_points_2d_p2
621 
mod_gauss_points_2d_p2
Definition: gauss_points_2d_p2.f90:4
def_type_mesh::mesh_type
Definition: def_type_mesh.f90:73
def_type_mesh
Definition: def_type_mesh.f90:4
mod_gauss_points_2d_p2::gauss_points_2d_p2
subroutine, public gauss_points_2d_p2(mesh)
Definition: gauss_points_2d_p2.f90:12
element_1d_p2
subroutine element_1d_p2(w, d, p)
Definition: gauss_points_2d_p2.f90:563
element_2d_p2_boundary
subroutine element_2d_p2_boundary(face, d, w)
Definition: gauss_points_2d_p2.f90:454
element_2d_p2
subroutine element_2d_p2(w, d, p)
Definition: gauss_points_2d_p2.f90:351