SFEMaNS  version 4.1 (work in progress)
Reference documentation for SFEMaNS
 All Classes Files Functions Variables Groups Pages
gauss_points_2d_p1.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_p1
5  PRIVATE
6 
7  PUBLIC gauss_points_2d_p1
8 
9 CONTAINS
10 
11  SUBROUTINE gauss_points_2d_p1(mesh)
12 
14  USE def_type_mesh
15 
16  IMPLICIT NONE
17 
18  TYPE(mesh_type), TARGET :: mesh ! > @param[out] 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 = 3, l_g = 3, &
26  n_ws = 2, l_gs = 2
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
49  REAL(KIND=8), DIMENSION(2) :: rnor, rsd ! JLG April 2009
50  !TEST
51  REAL(KIND=8), DIMENSION(n_w, l_G) :: ww_s
52  !TEST
53  !NEW Mai 26, 2004
54  INTEGER, DIMENSION(n_ws) :: nface
55  INTEGER, DIMENSION(l_Gs) :: ngauss
56  !NEW Mai 26, 2004
57  me => mesh%me
58  mes => mesh%mes
59  jj => mesh%jj
60  js => mesh%jjs
61  rr => mesh%rr
62 
63  ! JLG (April 2009)
64  IF (mesh%edge_stab) THEN
65  ALLOCATE(mesh%gauss%dwni(n_w, l_gs, 2, mesh%mi))
66  ALLOCATE(mesh%gauss%rji(l_gs, mesh%mi))
67  ALLOCATE(mesh%gauss%rnormsi(k_d, l_gs, 2, mesh%mi)) !JLG June 4 2012
68  ALLOCATE(mesh%gauss%wwsi(n_ws, l_gs)) !JLG June 4 2012
69  END IF
70  ! JLG (April 2009)
71 
72  ALLOCATE(mesh%gauss%ww(n_w, l_g))
73  ALLOCATE(mesh%gauss%wws(n_ws, l_gs))
74  ALLOCATE(mesh%gauss%dw(k_d, n_w, l_g, me ))
75  ALLOCATE(mesh%gauss%rnorms(k_d, l_gs, mes))
76  ALLOCATE(mesh%gauss%rj(l_g, me ))
77  ALLOCATE(mesh%gauss%rjs(l_gs, mes))
78  !NEW Juin 8, 2004
79  ALLOCATE(mesh%gauss%dw_s(k_d, n_w, l_gs, mesh%mes))
80  !NEW Juin 8, 2004
81  ALLOCATE(mesh%gauss%dwps(1, n_ws, l_gs, mes))
82  ALLOCATE(mesh%gauss%dws(1, n_ws, l_gs, mes))
83 
84  ww => mesh%gauss%ww
85  wws => mesh%gauss%wws
86  dw => mesh%gauss%dw
87  rnorms => mesh%gauss%rnorms
88  rj => mesh%gauss%rj
89  rjs => mesh%gauss%rjs
90  !NEW Mai 26, 2004
91  dw_s => mesh%gauss%dw_s
92  !NEW Mai 26, 2004
93  dwps => mesh%gauss%dwps
94  dws => mesh%gauss%dws
95 
96  mesh%gauss%k_d = k_d
97  mesh%gauss%n_w = n_w
98  mesh%gauss%l_G = l_g
99  mesh%gauss%n_ws = n_ws
100  mesh%gauss%l_Gs = l_gs
101 
102  ! evaluate and store the values of derivatives and of the
103  ! jacobian determinant at Gauss points of all volume elements
104 
105  ! volume elements
106 
107  CALL element_2d(ww, dd, pp)
108 
109  DO m = 1, me
110 
111  DO l = 1, l_g
112 
113  DO k = 1, k_d
114  r = rr(k, jj(:,m))
115  DO k1 = 1, k_d
116  dr(k, k1) = sum(r * dd(k1,:,l))
117  ENDDO
118  ENDDO
119 
120  rjac = dr(1,1)*dr(2,2) - dr(1,2)*dr(2,1)
121 
122  DO n = 1, n_w
123  dw(1, n, l, m) &
124  = (+ dd(1,n,l)*dr(2,2) - dd(2,n,l)*dr(2,1))/rjac
125  dw(2, n, l, m) &
126  = (- dd(1,n,l)*dr(1,2) + dd(2,n,l)*dr(1,1))/rjac
127  ENDDO
128 
129  rj(l, m) = rjac * pp(l)
130 
131  ENDDO
132 
133  ENDDO
134 
135  CALL element_1d(wws, dds, pps)
136 
137  ! surface elements
138  DO ms = 1, mes
139  m = mesh%neighs(ms)
140  !Il faut savoir quelle face est la bonne
141  DO n = 1, n_w
142  IF (minval(abs(js(:,ms)-jj(n,m)))==0) cycle
143  face = n
144  END DO
145  CALL element_2d_boundary(face, dd, ww_s)
146  n1 = modulo(face,3)+1; n2 = modulo(face+1,3)+1
147  IF (js(1,ms)==jj(n1,m) .AND. js(2,ms)==jj(n2,m)) THEN
148  nface(1)=n1; nface(2)=n2
149  ngauss(1)=1; ngauss(2)=2
150  ELSE IF (js(1,ms)==jj(n2,m) .AND. js(2,ms)==jj(n1,m)) THEN
151  nface(1)=n2; nface(2)=n1
152  ngauss(1)=2; ngauss(2)=1
153  ELSE
154  WRITE(*,*) ' BUG : gauss_points_2d_p1'
155  stop
156  END IF
157 
158  DO ls = 1, l_gs
159  !write(*,*) wws(:,ls)
160  !write(*,*) ww_s(nface,ngauss(ls))
161  DO k = 1, k_d
162  r = rr(k, jj(:,m))
163  DO k1 = 1, k_d
164  dr(k, k1) = sum(r * dd(k1,:,ls))
165  ENDDO
166  ENDDO
167 
168  rjac = (dr(1,1)*dr(2,2) - dr(1,2)*dr(2,1))
169 
170  DO n = 1, n_w
171  dw_s(1, n, ngauss(ls), ms) &
172  = (+ dd(1,n,ls)*dr(2,2) - dd(2,n,ls)*dr(2,1))/rjac
173  dw_s(2, n, ngauss(ls), ms) &
174  = (- dd(1,n,ls)*dr(1,2) + dd(2,n,ls)*dr(1,1))/rjac
175 
176  !dw_s(1, n, ls, ms) &
177  ! = (+ dd(1,nface(n),ngauss(ls))*dr(2,2) - dd(2,nface(n),ngauss(ls))*dr(2,1))/rjac
178  !dw_s(2, n, ls, ms) &
179  ! = (- dd(1,nface(n),ngauss(ls))*dr(1,2) + dd(2,nface(n),ngauss(ls))*dr(1,1))/rjac
180  !write(*,*) ms, n, dw_s(1, n, ls, ms), dw_s(2, n, ls, ms), mesh%sides(ms)
181  !write(*,*) ms, n, rr(1,js(n,ms)), rr(2,js(n,ms))
182 
183  ENDDO
184 
185  ENDDO
186 
187  ENDDO
188 
189  !CALL element_1d(wws, dds, pps)
190  DO ms = 1, mes
191 
192  DO ls = 1, l_gs
193 
194  DO k = 1, k_d
195  rs = rr(k, js(:,ms))
196  drs(1, k) = sum(rs * dds(1,:,ls))
197  ENDDO
198 
199  rjacs = sqrt( drs(1,1)**2 + drs(1,2)**2 )
200 
201  rnorms(1, ls, ms) = - drs(1,2)/rjacs ! outward normal
202  rnorms(2, ls, ms) = + drs(1,1)/rjacs ! outward normal
203  !erreur de signe corrigee le 23 mai 2001. La normale etait interieure
204  rjs(ls, ms) = rjacs * pps(ls)
205 
206  m = mesh%neighs(ms)
207  !Il faut savoir quelle face est la bonne
208  DO n = 1, n_w
209  IF (minval(abs(js(:,ms)-jj(n,m)))==0) cycle
210  face = n
211  END DO
212  rs = rr(:,jj(face,m)) - (rr(:,js(1,ms))+rr(:,js(2,ms)))/2
213  x = sum(rnorms(:,ls,ms)*rs)
214  IF (x>0) THEN
215  rnorms(:,ls,ms) = -rnorms(:,ls,ms)
216  !write(*,*) 'retournement de n', ms, ls
217  !erreur de signe corrigee le 7 juin 2004.
218  END IF
219 
220  ENDDO
221 
222  ENDDO
223 
224 
225  DO ms = 1, mes ! necessary only for evaluating gradient
226  ! tangential to the surface (ex COMMON Gauss_tan)
227 
228  DO ls = 1, l_gs
229  !TEST
230  !dwps(1, :, ls, ms) = dds(1, :, ls) / (rjs(ls, ms)/pps(ls))
231  !do n = 1, n_ws
232  !rjac = SUM(rnorms(:,ls,ms)*dw_s(:,n,ls,ms))
233  !write(*,*) n, ' js ', js(:,ms), rjs(ls, ms)/pps(ls)
234  !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)
235  !write(*,*) ms, ls, dwps(1, n, ls, ms)*rnorms(2,ls,ms), -dwps(1, n, ls, ms)*rnorms(1,ls,ms), &
236  ! SQRT(SUM(rnorms(:,ls,ms)**2)),rnorms(:,ls,ms)
237  !write(*,*) ' '
238  !end do
239  !TEST
240  dwps(1, :, ls, ms) = dds(1, :, ls) * pps(ls)
241  dws(1, :, ls, ms) = dds(1, :, ls)
242 
243  ENDDO
244 
245  ENDDO
246 
247  !-----------------------------------------------------------------------------
248  ! CELL INTERFACES (JLG, April 2009)
249  IF (mesh%edge_stab) THEN
250  CALL element_1d(mesh%gauss%wwsi, dds, pps) !JLG June 4 2012
251  DO ms = 1, mesh%mi
252  DO cote = 1, 2
253  m = mesh%neighi(cote,ms)
254  !Il faut savoir quelle face est la bonne
255  DO n = 1, n_w
256  IF (minval(abs(mesh%jjsi(:,ms)-jj(n,m)))==0) cycle
257  face = n
258  END DO
259  CALL element_2d_boundary(face, dd, ww_s)
260  n1 = modulo(face,3)+1; n2 = modulo(face+1,3)+1
261  IF (mesh%jjsi(1,ms)==jj(n1,m) .AND. mesh%jjsi(2,ms)==jj(n2,m)) THEN
262  nface(1)=n1; nface(2)=n2
263  ngauss(1)=1; ngauss(2)=2
264  ELSE IF (mesh%jjsi(1,ms)==jj(n2,m) .AND. mesh%jjsi(2,ms)==jj(n1,m)) THEN
265  nface(1)=n2; nface(2)=n1
266  ngauss(1)=2; ngauss(2)=1
267  ELSE
268  WRITE(*,*) ' BUG : gauss_points_2d_p1'
269  stop
270  END IF
271  DO ls = 1, l_gs
272  ! Compute normal vector
273  DO k = 1, k_d
274  rs = rr(k, mesh%jjsi(:,ms))
275  drs(1, k) = sum(rs * dds(1,:,ls))
276  ENDDO
277  rjacs = sqrt( drs(1,1)**2 + drs(1,2)**2 )
278  mesh%gauss%rji(ls, ms) = rjacs * pps(ls)
279  rnor(1) = drs(1,2)/rjacs
280  rnor(2) = - drs(1,1)/rjacs
281  mesh%gauss%rnormsi(1, ls, cote, ms) = rnor(1) !JLG, June 4 2012
282  mesh%gauss%rnormsi(2, ls, cote, ms) = rnor(2) !JLG, June 4 2012
283  rsd = rr(:,jj(face,m)) - (rr(:,mesh%jjsi(1,ms))+rr(:,mesh%jjsi(2,ms)))/2
284  x = sum(rnor(:)*rsd)
285  IF (x>0) THEN ! Verify the normal is outward
286  rnor = -rnor
287  !WRITE(*,*) 'retournement de n', ms, ls, cote, x
288  END IF
289  ! End computation normal vector
290  !WRITE(*,*) SUM((rr(:,mesh%jjsi(1,ms))-rr(:,mesh%jjsi(2,ms)))*rnor)
291  DO k = 1, k_d
292  r = rr(k, jj(:,m))
293  DO k1 = 1, k_d
294  dr(k, k1) = sum(r * dd(k1,:,ls))
295  ENDDO
296  ENDDO
297  rjac = (dr(1,1)*dr(2,2) - dr(1,2)*dr(2,1))
298  DO n = 1, n_w
299  !mesh%gauss%dwi(1, n, ngauss(ls), cote, ms) &
300  ! = (+ dd(1,n,ls)*dr(2,2) - dd(2,n,ls)*dr(2,1))/rjac
301  !mesh%gauss%dwi(2, n, ngauss(ls), cote, ms) &
302  ! = (- dd(1,n,ls)*dr(1,2) + dd(2,n,ls)*dr(1,1))/rjac
303  mesh%gauss%dwni(n, ngauss(ls), cote, ms) &
304  = rnor(1)*(+ dd(1,n,ls)*dr(2,2) - dd(2,n,ls)*dr(2,1))/rjac &
305  + rnor(2)*(- dd(1,n,ls)*dr(1,2) + dd(2,n,ls)*dr(1,1))/rjac
306  ENDDO
307  ENDDO
308  END DO
309  ENDDO
310 
311  DO ms = 1, mesh%mi
312  DO cote = 1, 2
313  m = mesh%neighi(cote,ms)
314  DO ls = 1, l_gs
315  x = abs(sqrt(sum(mesh%gauss%rnormsi(:, ls, cote, ms)**2))-1.d0)
316  !write(*,*) ' x ',x, ABS(SQRT(SUM(mesh%gauss%rnormsi(:, ls, cote, ms)**2)))
317  IF (x.GE.1.d-13) THEN
318  write(*,*) 'Bug1 in normal '
319  stop
320  END IF
321  rjac = (rr(1, mesh%jjsi(2,ms)) - rr(1, mesh%jjsi(1,ms)))*mesh%gauss%rnormsi(1, ls, cote, ms) &
322  + (rr(2, mesh%jjsi(2,ms)) - rr(2, mesh%jjsi(1,ms)))*mesh%gauss%rnormsi(2, ls, cote, ms)
323  x = sqrt((rr(1, mesh%jjsi(2,ms)) - rr(1, mesh%jjsi(1,ms)))**2 +&
324  (rr(2, mesh%jjsi(2,ms)) - rr(2, mesh%jjsi(1,ms)))**2)
325  !write(*,*) ' rjac/x ', rjac/x, x, rjac
326  IF (abs(rjac/x ).GE. 1.d-13) THEN
327  write(*,*) 'Bug2 in normal '
328  stop
329  END IF
330  END DO
331  END DO
332  END DO
333  END IF
334  !-----------------------------------------------------------------------------
335 
336  !TEST
337  !WRITE(*,*) MAXVAL(ABS(dw(:,:,1,mesh%neighi(1,:)) - mesh%gauss%dwi(:,:,1,1,:)))
338  !WRITE(*,*) MAXVAL(ABS(dw(:,:,2,mesh%neighi(2,:)) - mesh%gauss%dwi(:,:,1,2,:)))
339  !WRITE(*,*) MAXVAL(ABS(jj(:,mesh%neighi(1,:))-mesh%jji(:,1,:)))
340  !WRITE(*,*) MAXVAL(ABS(jj(:,mesh%neighi(2,:))-mesh%jji(:,2,:)))
341  !TEST
342 
343  !PRINT*, 'end of gen_Gauss'
344 
345  CONTAINS
346 
347  SUBROUTINE element_2d (w, d, p)
348 
352 
353  ! w(n_w, l_G) : values of shape functions at Gauss points
354  ! d(2, n_w, l_G) : derivatives values of shape functions at Gauss points
355  ! p(l_G) : weight for Gaussian quadrature at Gauss points
356 
357  IMPLICIT NONE
358 
359  REAL(KIND=8), DIMENSION( n_w, l_G), INTENT(OUT) :: w
360  REAL(KIND=8), DIMENSION(2, n_w, l_G), INTENT(OUT) :: d
361  REAL(KIND=8), DIMENSION(l_G), INTENT(OUT) :: p
362 
363  REAL(KIND=8), DIMENSION(l_G) :: xx, yy
364  INTEGER :: j
365 
366  REAL(KIND=8) :: zero = 0, one = 1, two = 2, three = 3, six = 6
367 
368  REAL(KIND=8) :: f1, f2, f3, x, y
369  f1(x, y) = one - x - y
370  f2(x, y) = x
371  f3(x, y) = y
372 
373  xx(1) = one/six; xx(2) = two/three; xx(3) = one/six
374  yy(1) = one/six; yy(2) = one/six; yy(3) = two/three
375 
376  DO j = 1, l_g
377 
378  w(1, j) = f1(xx(j), yy(j))
379  d(1, 1, j) = - one
380  d(2, 1, j) = - one
381 
382  w(2, j) = f2(xx(j), yy(j))
383  d(1, 2, j) = one
384  d(2, 2, j) = zero
385 
386  w(3, j) = f3(xx(j), yy(j))
387  d(1, 3, j) = zero
388  d(2, 3, j) = one
389 
390  p(j) = one/six
391 
392  ENDDO
393 
394  END SUBROUTINE element_2d
395 
396  SUBROUTINE element_2d_boundary (face, d, w)
397 
401 
402  ! face : face
403  ! d(2, n_w, l_G) : derivatives values of shape functions at Gauss points
404  ! w(n_w, l_G) : values of shape functions at Gauss points
405 
406 
407 
408  IMPLICIT NONE
409 
410  INTEGER, INTENT(IN) :: face
411  REAL(KIND=8), DIMENSION(2, n_w, l_Gs), INTENT(OUT) :: d
412  REAL(KIND=8), DIMENSION( n_w, l_G), INTENT(OUT) :: w
413 
414  REAL(KIND=8), DIMENSION(l_Gs) :: xx, yy
415  INTEGER :: j
416 
417  REAL(KIND=8) :: zero = 0, one = 1, three=3, half = 0.5d0
418 
419  REAL(KIND=8) :: f1, f2, f3, x, y
420  f1(x, y) = 1-x-y
421  f2(x, y) = x
422  f3(x, y) = y
423 
424  IF (face==1) THEN
425  xx(1) = half + half/sqrt(three)
426  yy(1) = half - half/sqrt(three)
427 
428  xx(2) = half - half/sqrt(three)
429  yy(2) = half + half/sqrt(three)
430  ELSE IF (face==2) THEN
431  xx(1) = zero
432  yy(1) = half + half/sqrt(three)
433 
434  xx(2) = zero
435  yy(2) = half - half/sqrt(three)
436  ELSE
437  xx(1) = half - half/sqrt(three)
438  yy(1) = zero
439 
440  xx(2) = half + half/sqrt(three)
441  yy(2) = zero
442  END IF
443 
444  DO j = 1, l_gs
445  w(1,j) = f1(xx(j), yy(j))
446  d(1, 1, j) = -one
447  d(2, 1, j) = -one
448 
449  w(2,j) = f2(xx(j), yy(j))
450  d(1, 2, j) = one
451  d(2, 2, j) = zero
452 
453  w(3,j) = f3(xx(j), yy(j))
454  d(1, 3, j) = zero
455  d(2, 3, j) = one
456  ENDDO
457 
458  END SUBROUTINE element_2d_boundary
459 
460  !------------------------------------------------------------------------------
461 
462  SUBROUTINE element_1d (w, d, p)
463 
464  ! one-dimensional element with linear interpolation
465  ! and two Gauss integration points
466 
467  ! w(n_w, l_G) : values of shape functions at Gauss points
468  ! d(1, n_w, l_G) : derivatives values of shape functions at Gauss points
469  ! p(l_G) : weight for Gaussian quadrature at Gauss points
470 
471  IMPLICIT NONE
472 
473  REAL(KIND=8), DIMENSION( n_ws, l_Gs), INTENT(OUT) :: w
474  REAL(KIND=8), DIMENSION(1, n_ws, l_Gs), INTENT(OUT) :: d
475  REAL(KIND=8), DIMENSION(l_Gs), INTENT(OUT) :: p
476 
477  REAL(KIND=8), DIMENSION(l_Gs) :: xx
478  INTEGER :: j
479 
480  REAL(KIND=8) :: one = 1, two = 2, three = 3
481 
482  REAL(KIND=8) :: f1, f2, x
483  f1(x) = (one - x)/two
484  f2(x) = (x + one)/two
485 
486  xx(1) = - one/sqrt(three)
487  xx(2) = + one/sqrt(three)
488 
489  DO j = 1, l_gs
490 
491  w(1, j) = f1(xx(j))
492  d(1, 1, j) = - one/two
493 
494  w(2, j) = f2(xx(j))
495  d(1, 2, j) = + one/two
496 
497  p(j) = one
498 
499  ENDDO
500 
501  END SUBROUTINE element_1d
502 
503 
504  END SUBROUTINE gauss_points_2d_p1
505 
506 END MODULE mod_gauss_points_2d_p1
507 
element_1d
subroutine element_1d(w, d, p)
Definition: gauss_points_2d_p1.f90:462
def_type_mesh::mesh_type
Definition: def_type_mesh.f90:73
element_2d_boundary
subroutine element_2d_boundary(face, d, w)
Defines gauss quadrature for faces of triangular element with linear interpolation and two Gauss inte...
Definition: gauss_points_2d_p1.f90:396
element_2d
subroutine element_2d(w, d, p)
Defines gauss quadrature for triangular element with linear interpolation and three Gauss integration...
Definition: gauss_points_2d_p1.f90:347
mod_gauss_points_2d_p1
Definition: gauss_points_2d_p1.f90:4
def_type_mesh
Definition: def_type_mesh.f90:4
mod_gauss_points_2d_p1::gauss_points_2d_p1
subroutine, public gauss_points_2d_p1(mesh)
Prepares gauss points for P1 approximation in 2D.
Definition: gauss_points_2d_p1.f90:11