SFEMaNS  version 5.3
Reference documentation for SFEMaNS
gauss_points_2d.f90
Go to the documentation of this file.
1 !===
2 !Author: Jean-Luc Guermond, Copyright December 24th 2018
3 !===
4 MODULE mod_gauss_points_2d
5  PRIVATE
6  PUBLIC gauss_points_2d
7  INTEGER, PARAMETER :: k_d=2
8 CONTAINS
9  SUBROUTINE gauss_points_2d(mesh,type_fe)
10  USE def_type_mesh
11  USE gp_2d_p1
12  USE gp_2d_p2
13  USE gp_2d_p3
14  USE sub_plot
15  IMPLICIT NONE
16  TYPE(mesh_type), TARGET :: mesh
17  INTEGER :: type_fe
18  INTEGER, POINTER :: me, mes
19  INTEGER, DIMENSION(:,:), POINTER :: jj
20  INTEGER, DIMENSION(:,:), POINTER :: js
21  REAL(KIND=8), DIMENSION(:,:), POINTER :: rr
22  REAL(KIND=8), DIMENSION(:,:), POINTER :: ww
23  REAL(KIND=8), DIMENSION(:,:), POINTER :: wws
24  REAL(KIND=8), DIMENSION(:,:,:,:), POINTER :: dw
25  REAL(KIND=8), DIMENSION(:,:,:), POINTER :: rnorms
26  REAL(KIND=8), DIMENSION(:,:), POINTER :: rj
27  REAL(KIND=8), DIMENSION(:,:), POINTER :: rjs
28  REAL(KIND=8), DIMENSION(:,:,:,:), POINTER :: dw_s
29  REAL(KIND=8), DIMENSION(:,:,:,:), POINTER :: dwps
30  REAL(KIND=8), DIMENSION(:,:,:,:), POINTER :: dws
31  REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: dd
32  REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: dds
33  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dds_v
34  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: pp
35  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: pps
36  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: r
37  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: rs
38  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: ww_s
39  INTEGER, DIMENSION(:), ALLOCATABLE :: ngauss
40  REAL(KIND=8), DIMENSION(k_d, k_d) :: dr
41  REAL(KIND=8), DIMENSION(1, k_d) :: drs
42  REAL(KIND=8) :: rjac, rjacs, x
43  INTEGER :: m, l, k, k1, n, n1, n2, ms, ns, ls, face, cote, orient
44  INTEGER :: n_w, n_ws, l_G, l_Gs
45  REAL(KIND=8), DIMENSION(k_d) :: rnor, rsd
46 
47  SELECT CASE(type_fe)
48  CASE(1)
49  n_w=3; n_ws=2; l_g=3; l_gs=2
50  CASE(2)
51  n_w=6; n_ws=3; l_g=7; l_gs=3
52  CASE(3)
53  n_w=10; n_ws=4; l_g=7; l_gs=4
54  CASE DEFAULT
55  WRITE(*,*) ' FE not programmed yet', type_fe
56  stop
57  END SELECT
58  ALLOCATE(dd(k_d,n_w,l_g),dds(1,n_ws,l_gs),pp(l_g),pps(l_gs),&
59  r(n_w),rs(n_ws),ww_s(n_w,l_gs),dds_v(n_ws,n_ws),ngauss(l_gs))
60  me => mesh%me
61  mes => mesh%mes
62  jj => mesh%jj
63  js => mesh%jjs
64  rr => mesh%rr
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))
69  ALLOCATE(mesh%gauss%wwsi(n_ws, l_gs))
70  END IF
71  ALLOCATE(mesh%gauss%ww(n_w, l_g))
72  ALLOCATE(mesh%gauss%wws(n_ws, l_gs))
73  ALLOCATE(mesh%gauss%dw(k_d, n_w, l_g, me ))
74  ALLOCATE(mesh%gauss%rnorms(k_d, l_gs, mes))
75  ALLOCATE(mesh%gauss%rnorms_v(k_d,n_ws, mes))
76  ALLOCATE(mesh%gauss%rj(l_g, me ))
77  ALLOCATE(mesh%gauss%rjs(l_gs, mes))
78  ALLOCATE(mesh%gauss%dw_s(k_d, n_w, l_gs, mesh%mes))
79  ALLOCATE(mesh%gauss%dwps(1, n_ws, l_gs, mes))
80  ALLOCATE(mesh%gauss%dws(1, n_ws, l_gs, mes))
81  ww => mesh%gauss%ww
82  wws => mesh%gauss%wws
83  dw => mesh%gauss%dw
84  rnorms => mesh%gauss%rnorms
85  rj => mesh%gauss%rj
86  rjs => mesh%gauss%rjs
87  dw_s => mesh%gauss%dw_s
88  dwps => mesh%gauss%dwps
89  dws => mesh%gauss%dws
90  mesh%gauss%k_d = k_d
91  mesh%gauss%n_w = n_w
92  mesh%gauss%l_G = l_g
93  mesh%gauss%n_ws = n_ws
94  mesh%gauss%l_Gs = l_gs
95  !===Evaluate and store the values of derivatives and of the
96  !===Jacobian determinant at Gauss points of all volume elements
97  SELECT CASE(type_fe)
98  CASE(1)
99  CALL element_2d_p1(ww, dd, pp, n_w, l_g)
100  CASE(2)
101  CALL element_2d_p2(ww, dd, pp, n_w, l_g)
102  CASE(3)
103  CALL element_2d_p3(ww, dd, pp, n_w, l_g)
104  END SELECT
105  DO m = 1, me
106  DO l = 1, l_g
107  DO k = 1, k_d
108  r = rr(k, jj(:,m))
109  DO k1 = 1, k_d
110  dr(k, k1) = sum(r * dd(k1,:,l))
111  ENDDO
112  ENDDO
113  rjac = dr(1,1)*dr(2,2) - dr(1,2)*dr(2,1)
114  DO n = 1, n_w
115  dw(1, n, l, m) &
116  = (+ dd(1,n,l)*dr(2,2) - dd(2,n,l)*dr(2,1))/rjac
117  dw(2, n, l, m) &
118  = (- dd(1,n,l)*dr(1,2) + dd(2,n,l)*dr(1,1))/rjac
119  ENDDO
120  rj(l, m) = abs(rjac)*pp(l) !===Sign of Jacobian determinant unknown
121  ENDDO
122  ENDDO
123  SELECT CASE(type_fe)
124  CASE(1)
125  CALL element_1d_p1(wws, dds, pps, n_ws, l_gs)
126  CASE(2)
127  CALL element_1d_p2(wws, dds, pps, n_ws, l_gs)
128  CASE(3)
129  CALL element_1d_p3(wws, dds, pps, n_ws, l_gs)
130  END SELECT
131  !===surface elements
132  DO ms = 1, mes
133  m = mesh%neighs(ms)
134  !===Determine which face of the reference element is associated with ms
135  DO n = 1, 3
136  IF (minval(abs(js(:,ms)-jj(n,m)))==0) cycle
137  face = n
138  END DO
139  SELECT CASE(type_fe)
140  CASE(1)
141  CALL element_2d_p1_boundary(face, dd, ww_s, n_w, l_gs)
142  CASE(2)
143  CALL element_2d_p2_boundary(face, dd, ww_s, n_w, l_gs)
144  CASE(3)
145  CALL element_2d_p3_boundary(face, dd, ww_s, n_w, l_gs)
146  END SELECT
147  !===Fix for meshes not oriented correctly
148  n1 = modulo(face,3)+1; n2 = modulo(face+1,3)+1
149  IF (js(1,ms)==jj(n1,m) .AND. js(2,ms)==jj(n2,m)) THEN
150  orient = 1
151  ELSE
152  orient = -1
153  END IF
154  IF (face==2) THEN
155  orient = -orient
156  END IF
157  !=========TEST
158  !orient = 1
159  !=========TEST
160  DO ls = 1, l_gs
161  ngauss(ls) = orient*ls+(l_gs+1)*(1-orient)/2 !===rearrangement of gauss points
162  END DO
163  !===End fix
164  DO ls = 1, l_gs
165  DO k = 1, k_d
166  r = rr(k, jj(:,m))
167  DO k1 = 1, k_d
168  dr(k, k1) = sum(r * dd(k1,:,ls))
169  ENDDO
170  ENDDO
171  rjac = (dr(1,1)*dr(2,2) - dr(1,2)*dr(2,1))
172  DO n = 1, n_w
173  dw_s(1,n,ngauss(ls),ms) = (+dd(1,n,ls)*dr(2,2) - dd(2,n,ls)*dr(2,1))/rjac
174  dw_s(2,n,ngauss(ls),ms) = (-dd(1,n,ls)*dr(1,2) + dd(2,n,ls)*dr(1,1))/rjac
175  ENDDO
176  ENDDO
177  ENDDO
178  DO ms = 1, mes
179  DO ls = 1, l_gs
180  DO k = 1, k_d
181  rs = rr(k, js(:,ms))
182  drs(1, k) = sum(rs*dds(1,:,ls))
183  ENDDO
184  rjacs = sqrt(drs(1,1)**2 + drs(1,2)**2)
185  rnorms(1, ls, ms) = - drs(1,2)/rjacs
186  rnorms(2, ls, ms) = + drs(1,1)/rjacs
187  rjs(ls, ms) = rjacs * pps(ls)
188  m = mesh%neighs(ms)
189  DO n = 1, n_w
190  IF (minval(abs(js(:,ms)-jj(n,m)))==0) cycle
191  face = n
192  END DO
193  rsd(1:k_d) = rr(:,jj(face,m)) - (rr(:,js(1,ms))+rr(:,js(2,ms)))/2
194  x = sum(rnorms(:,ls,ms)*rsd(1:k_d))
195  IF (x>0) THEN
196  rnorms(:,ls,ms) = -rnorms(:,ls,ms)
197  END IF
198  ENDDO
199  ENDDO
200  DO ms = 1, mes !===Evaluation of tangent gradient. Obsolete
201  DO ls = 1, l_gs
202  dwps(1,:,ls,ms) = dds(1,:,ls)*pps(ls)
203  dws(1,:,ls,ms) = dds(1,:,ls)
204  ENDDO
205  ENDDO
206 
207  !===Normal vector at Lagrange nodes
208  SELECT CASE(type_fe)
209  CASE(1)
210  CALL element_1d_p1_at_nodes (dds_v, n_ws)
211  CASE DEFAULT
212  !WRITE(*,*) 'BUG gauss_points_2d, rnorms_v not programmed yet'
213  END SELECT
214  DO ms = 1, mes
215  DO ns = 1, n_ws
216  DO k = 1, k_d
217  rs = rr(k, js(:,ms))
218  drs(1, k) = sum(rs * dds_v(:,ns))
219  ENDDO
220  rjacs = sqrt( drs(1,1)**2 + drs(1,2)**2 )
221  mesh%gauss%rnorms_v(1, ns, ms) = -drs(1,2)/rjacs
222  mesh%gauss%rnorms_v(2, ns, ms) = drs(1,1)/rjacs
223  m = mesh%neighs(ms)
224  !===Find correct face and orient normal outwards
225  DO n = 1, n_w
226  IF (minval(abs(js(:,ms)-jj(n,m)))==0) cycle
227  face = n
228  END DO
229  rs = rr(:,jj(face,m)) - (rr(:,js(1,ms))+rr(:,js(2,ms)))/2
230  x = sum(mesh%gauss%rnorms_v(:,ns,ms)*rs)
231  IF (x>0) THEN
232  mesh%gauss%rnorms_v(:,ns,ms) = - mesh%gauss%rnorms_v(:,ns,ms)
233  END IF
234  END DO
235  ENDDO
236 
237  !===Cell interface (JLG, April 2009)
238  IF (mesh%edge_stab) THEN
239  SELECT CASE(type_fe)
240  CASE(1)
241  CALL element_1d_p1(mesh%gauss%wwsi, dds, pps, n_ws, l_gs)
242  CASE(2)
243  CALL element_1d_p2(mesh%gauss%wwsi, dds, pps, n_ws, l_gs)
244  CASE(3)
245  CALL element_1d_p3(mesh%gauss%wwsi, dds, pps, n_ws, l_gs)
246  END SELECT
247  DO ms = 1, mesh%mi
248  DO cote = 1, 2
249  m = mesh%neighi(cote,ms)
250  DO n = 1, 3
251  IF (minval(abs(mesh%jjsi(1:2,ms)-jj(n,m)))==0) cycle
252  face = n
253  END DO
254  SELECT CASE(type_fe)
255  CASE(1)
256  CALL element_2d_p1_boundary(face, dd, ww_s, n_w, l_gs)
257  CASE(2)
258  CALL element_2d_p2_boundary(face, dd, ww_s, n_w, l_gs)
259  CASE(3)
260  CALL element_2d_p3_boundary(face, dd, ww_s, n_w, l_gs)
261  END SELECT
262  !===Fix for meshes not oriented correctly
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  orient = 1
266  ELSE
267  orient = -1
268  END IF
269  IF (face==2) THEN
270  orient = -orient
271  END IF
272  DO ls = 1, l_gs
273  ngauss(ls) = orient*ls+(l_gs+1)*(1-orient)/2 !===rearrangement of gauss points
274  END DO
275  !===End fix
276  DO ls = 1, l_gs
277  DO k = 1, k_d
278  rs = rr(k,mesh%jjsi(:,ms))
279  drs(1, k) = sum(rs*dds(1,:,ls))
280  ENDDO
281  rjacs = sqrt(drs(1,1)**2 + drs(1,2)**2)
282  mesh%gauss%rji(ls, ms) = rjacs*pps(ls)
283  rnor(1) = drs(1,2)/rjacs
284  rnor(2) = - drs(1,1)/rjacs
285  mesh%gauss%rnormsi(1, ls, cote, ms) = rnor(1)
286  mesh%gauss%rnormsi(2, ls, cote, ms) = rnor(2)
287  rsd = rr(:,jj(face,m)) - (rr(:,mesh%jjsi(1,ms))+rr(:,mesh%jjsi(2,ms)))/2
288  x = sum(rnor(:)*rsd)
289  IF (x>0) THEN !===Verify the normal is outward
290  rnor = -rnor
291  END IF
292  DO k = 1, k_d
293  r = rr(k, jj(:,m))
294  DO k1 = 1, k_d
295  dr(k, k1) = sum(r * dd(k1,:,ls))
296  ENDDO
297  ENDDO
298  rjac = (dr(1,1)*dr(2,2) - dr(1,2)*dr(2,1))
299  DO n = 1, n_w !===Fix with ngauss
300  mesh%gauss%dwni(n, ngauss(ls), cote, ms) &
301  = rnor(1)*(+ dd(1,n,ls)*dr(2,2) - dd(2,n,ls)*dr(2,1))/rjac &
302  + rnor(2)*(- dd(1,n,ls)*dr(1,2) + dd(2,n,ls)*dr(1,1))/rjac
303  ENDDO
304  ENDDO
305  END DO
306  ENDDO
307  !===Some verifications
308  DO ms = 1, mesh%mi
309  DO cote = 1, 2
310  m = mesh%neighi(cote,ms)
311  DO ls = 1, l_gs
312  x = abs(sqrt(sum(mesh%gauss%rnormsi(:, ls, cote, ms)**2))-1.d0)
313  IF (x.GE.1.d-13) THEN
314  write(*,*) 'Bug1 in normal '
315  stop
316  END IF
317  rjac = (rr(1, mesh%jjsi(2,ms))-rr(1, mesh%jjsi(1,ms)))*mesh%gauss%rnormsi(1, ls, cote, ms) &
318  + (rr(2, mesh%jjsi(2,ms))-rr(2, mesh%jjsi(1,ms)))*mesh%gauss%rnormsi(2, ls, cote, ms)
319  x = sqrt((rr(1, mesh%jjsi(2,ms))-rr(1, mesh%jjsi(1,ms)))**2 +&
320  (rr(2, mesh%jjsi(2,ms))-rr(2, mesh%jjsi(1,ms)))**2)
321  IF (abs(rjac/x) .GE. 1.d-10) THEN
322  write(*,*) 'Bug2 in normal ', ms, cote, ls, abs(rjac/x)
323  stop
324  END IF
325  END DO
326  END DO
327  END DO
328  END IF
329 
330  DEALLOCATE(dd,dds,pp,pps,r,rs,ww_s,dds_v,ngauss)
331  END SUBROUTINE gauss_points_2d
332 END MODULE mod_gauss_points_2d
333 
gp_2d_p1
Definition: GP_2D_p1.f90:5
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
sub_plot
Definition: sub_plot.f90:4
gp_2d_p1::element_1d_p1
subroutine, public element_1d_p1(w, d, p, n_ws, l_Gs)
Definition: GP_2D_p1.f90:93
gp_2d_p2::element_2d_p2
subroutine, public element_2d_p2(w, d, p, n_w, l_G)
Definition: GP_2D_p2.f90:10
mod_gauss_points_2d::k_d
integer, parameter k_d
Definition: gauss_points_2d.f90:7
mod_gauss_points_2d
Definition: gauss_points_2d.f90:4
gp_2d_p2
Definition: GP_2D_p2.f90:5
def_type_mesh::mesh_type
Definition: def_type_mesh.f90:74
gp_2d_p2::element_2d_p2_boundary
subroutine, public element_2d_p2_boundary(face, d, w, n_w, l_Gs)
Definition: GP_2D_p2.f90:91
gp_2d_p1::element_2d_p1
subroutine, public element_2d_p1(w, d, p, n_w, l_G)
Definition: GP_2D_p1.f90:10
gp_2d_p3
Definition: GP_2D_p3.f90:4
gp_2d_p2::element_1d_p2
subroutine, public element_1d_p2(w, d, p, n_ws, l_Gs)
Definition: GP_2D_p2.f90:177
gp_2d_p1::element_1d_p1_at_nodes
subroutine, public element_1d_p1_at_nodes(d, n_ws)
Definition: GP_2D_p1.f90:121
gp_2d_p1::element_2d_p1_boundary
subroutine, public element_2d_p1_boundary(face, d, w, n_w, l_Gs)
Definition: GP_2D_p1.f90:44
def_type_mesh
Definition: def_type_mesh.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
mod_gauss_points_2d::gauss_points_2d
subroutine, public gauss_points_2d(mesh, type_fe)
Definition: gauss_points_2d.f90:10