SFEMaNS  version 5.3
Reference documentation for SFEMaNS
fem_tn_axi.f90
Go to the documentation of this file.
1 !
2 !Authors: Jean-Luc Guermond, Raphael Laguerre, Luigi Quartapelle, Copyrights 1996, 2000, 2004
3 !
4 MODULE fem_tn_axi
5 CONTAINS
6 
7 
8 SUBROUTINE dot_product (mesh, ff, gg, t)
9 !===============================
10 
11 ! sqrt(< f^2 >) ===> t
12 
13  USE gauss_points
14 
15  IMPLICIT NONE
16 
17  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff, gg
18  REAL(KIND=8), INTENT(OUT) :: t
19 
20  INTEGER :: m, l ,i ,ni
21 
22  TYPE(mesh_type), TARGET :: mesh
23  INTEGER, DIMENSION(:,:), POINTER :: jj
24  INTEGER, POINTER :: me
25  REAL(KIND=8) ,DIMENSION(:,:), POINTER :: rr
26  REAL(KIND=8) :: ray
27 
28  CALL gauss(mesh)
29  jj => mesh%jj
30  me => mesh%me
31  rr=> mesh%rr
32  t = 0
33 
34  DO m = 1, me
35  DO l = 1, l_g
36 
37 !===Compute radius of Gauss point
38  ray = 0
39  DO ni = 1, n_w; i = jj(ni,m)
40  ray = ray + rr(1,i)*ww(ni,l)
41  END DO
42 
43  t = t + sum(ff(jj(:,m)) * ww(:,l))*sum(gg(jj(:,m)) * ww(:,l))*rj(l,m)*ray
44 
45  ENDDO
46  ENDDO
47 
48 END SUBROUTINE dot_product
49 
50 
51 SUBROUTINE ns_0 (mesh, ff, t)
52 !===============================
53 
54 ! sqrt(< f^2 >) ===> t
55 
56  USE gauss_points
57 
58  IMPLICIT NONE
59 
60  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff
61  REAL(KIND=8), INTENT(OUT) :: t
62 
63  INTEGER :: m, l ,i ,ni
64 
65  TYPE(mesh_type), TARGET :: mesh
66  INTEGER, DIMENSION(:,:), POINTER :: jj
67  INTEGER, POINTER :: me
68  REAL(KIND=8) ,DIMENSION(:,:), POINTER :: rr
69  REAL(KIND=8) :: ray
70 
71  CALL gauss(mesh)
72  jj => mesh%jj
73  me => mesh%me
74  rr=> mesh%rr
75  t = 0
76 
77  DO m = 1, me
78  DO l = 1, l_g
79 
80 !===Compute radius of Gauss point
81  ray = 0
82  DO ni = 1, n_w; i = jj(ni,m)
83  ray = ray + rr(1,i)*ww(ni,l)
84  END DO
85 
86  t = t + sum(ff(jj(:,m)) * ww(:,l))**2 * rj(l,m)*ray
87 
88  ENDDO
89  ENDDO
90 
91  t = sqrt(t)
92 
93 END SUBROUTINE ns_0
94 
95 SUBROUTINE ns_l1 (mesh, ff, t)
96 !===============================
97 
98 ! < f > ===> t
99 
100  USE gauss_points
101 
102  IMPLICIT NONE
103 
104  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff
105  REAL(KIND=8), INTENT(OUT) :: t
106 
107  INTEGER :: m, l ,i ,ni
108 
109  TYPE(mesh_type), TARGET :: mesh
110  INTEGER, DIMENSION(:,:), POINTER :: jj
111  INTEGER, POINTER :: me
112  REAL(KIND=8) ,DIMENSION(:,:), POINTER :: rr
113  REAL(KIND=8) :: ray
114 
115  CALL gauss(mesh)
116  jj => mesh%jj
117  me => mesh%me
118  rr=> mesh%rr
119  t = 0
120 
121  DO m = 1, me
122  DO l = 1, l_g
123 
124 !===Compute radius of Gauss point
125  ray = 0
126  DO ni = 1, n_w; i = jj(ni,m)
127  ray = ray + rr(1,i)*ww(ni,l)
128  END DO
129 
130  t = t + sum(ff(jj(:,m)) * ww(:,l))* rj(l,m)*ray
131 
132  ENDDO
133  ENDDO
134 
135 END SUBROUTINE ns_l1
136 
137 SUBROUTINE average (mesh, ff, t)
138  !===============================
139 
140  ! < f >/vol ===> t
141 
142  USE gauss_points
143 
144  IMPLICIT NONE
145 
146  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff
147  REAL(KIND=8), INTENT(OUT) :: t
148 
149  INTEGER :: m, l ,i ,ni
150 
151  TYPE(mesh_type), TARGET :: mesh
152  INTEGER, DIMENSION(:,:), POINTER :: jj
153  INTEGER, POINTER :: me
154  REAL(KIND=8) ,DIMENSION(:,:), POINTER :: rr
155  REAL(KIND=8) :: ray, vol
156 
157 
158  CALL gauss(mesh)
159  jj => mesh%jj
160  me => mesh%me
161  rr=> mesh%rr
162  t = 0
163  vol = 0.d0
164  IF (SIZE(ff)==mesh%np) THEN
165  DO m = 1, me
166  DO l = 1, l_g
167 
168  !===Compute radius of Gauss point
169  ray = 0
170  DO ni = 1, n_w; i = jj(ni,m)
171  ray = ray + rr(1,i)*ww(ni,l)
172  END DO
173 
174  t = t + sum(ff(jj(:,m)) * ww(:,l))* rj(l,m)*ray
175  vol = vol + rj(l,m)*ray
176  ENDDO
177  ENDDO
178  ELSE IF (SIZE(ff)==mesh%me) THEN
179  DO m = 1, me
180  DO l = 1, l_g
181 
182  !===Compute radius of Gauss point
183  ray = 0
184  DO ni = 1, n_w; i = jj(ni,m)
185  ray = ray + rr(1,i)*ww(ni,l)
186  END DO
187 
188  t = t + ff(m)* rj(l,m)*ray
189  vol = vol + rj(l,m)*ray
190  ENDDO
191  ENDDO
192  ELSE
193  WRITE(*,*) ' BUG in average '
194  stop
195  END IF
196  t = t / vol
197 
198 END SUBROUTINE average
199 
200 SUBROUTINE ns_anal_0 (mesh, ff, ff_anal, t)
201 !===============================
202 
203 ! sqrt(< f^2 >) ===> t
204 
205  USE gauss_points
206 
207  IMPLICIT NONE
208 
209  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff
210  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff_anal
211  REAL(KIND=8), INTENT(OUT) :: t
212 
213  INTEGER :: m, l ,i , ni, index
214 
215  TYPE(mesh_type), TARGET :: mesh
216  INTEGER, DIMENSION(:,:), POINTER :: jj
217  INTEGER, POINTER :: me
218  REAL(KIND=8) ,DIMENSION(:,:), POINTER :: rr
219  REAL(KIND=8) :: ray, fl
220 
221  CALL gauss(mesh)
222  jj => mesh%jj
223  me => mesh%me
224  rr=> mesh%rr
225  t = 0
226 
227  DO m = 1, me
228  index = (m-1)*l_g
229  DO l = 1, l_g; index = index + 1
230 
231 !===Compute radius of Gauss point
232  ray = 0
233  fl = 0
234  DO ni = 1, n_w; i = jj(ni,m)
235  ray = ray + rr(1,i)*ww(ni,l)
236  fl = fl + ff(i) * ww(ni,l)
237  END DO
238 
239  t = t + (fl - ff_anal(index))**2 * ray* rj(l,m)
240 
241  ENDDO
242  ENDDO
243 
244  t = sqrt(t)
245 
246 END SUBROUTINE ns_anal_0
247 
248 SUBROUTINE ns_1 (mesh, ff, t)
249 !===============================
250 
251 ! sqrt(<< (Df).(Df) >>) ===> t
252 
253  USE gauss_points
254 
255  IMPLICIT NONE
256 
257  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff
258  REAL(KIND=8), INTENT(OUT) :: t
259 
260  INTEGER :: m, l, k,i,ni
261  REAL(KIND=8) :: s
262 
263  TYPE(mesh_type), TARGET :: mesh
264  INTEGER, DIMENSION(:,:), POINTER :: jj
265  INTEGER, POINTER :: me
266  REAL(KIND=8) ,DIMENSION(:,:), POINTER :: rr
267  REAL(KIND=8) :: ray
268 
269  CALL gauss(mesh)
270  jj => mesh%jj
271  me => mesh%me
272  rr=> mesh%rr
273 
274  t = 0
275 
276  DO m = 1, me
277  DO l = 1, l_g
278  s = 0
279  DO k = 1, k_d
280 
281  !===Compute radius of Gauss point
282  ray = 0
283  DO ni = 1, n_w; i = jj(ni,m)
284  ray = ray + rr(1,i)*ww(ni,l)
285  END DO
286 
287  s = s + sum(ff(jj(:,m)) * dw(k,:,l,m))**2*ray
288 
289  ENDDO
290 
291  t = t + s * rj(l,m)
292 
293  ENDDO
294  ENDDO
295 
296  t = sqrt(t)
297 
298 END SUBROUTINE ns_1
299 
300 SUBROUTINE nv_1(mesh, mod_max, ff, t)
301 !semi-norme vectoriel H1 pour une representation de Fourier
302 !(v_rc, v_rs, v_thc, v_ths, v_zc, v_zs)
303 !sqrt(<< (Dv).(Dv) >>) ===> t
304 
305  USE gauss_points
306 
307  IMPLICIT NONE
308 
309  TYPE(mesh_type), TARGET :: mesh
310  INTEGER, INTENT(IN) :: mod_max
311  REAL(KIND=8), DIMENSION(6,mesh%np,0:mod_max), INTENT(IN) :: ff
312  REAL(KIND=8), INTENT(OUT) :: t
313 
314  INTEGER :: m, l,i,ni ,j
315 
316  INTEGER, DIMENSION(:,:), POINTER :: jj
317  INTEGER, POINTER :: me
318  REAL(KIND=8) ,DIMENSION(:,:), POINTER :: rr
319  REAL(KIND=8) :: ray
320  REAL(KIND=8), DIMENSION(2) :: div
321 
322  CALL gauss(mesh)
323  jj => mesh%jj
324  me => mesh%me
325  rr=> mesh%rr
326 
327  t = 0.d0
328  div = 0.d0
329 
330  DO j = 0, mod_max
331  DO m = 1, me
332  DO l = 1, l_g
333 !===Compute radius of Gauss point
334  ray = 0
335  DO ni = 1, n_w; i = jj(ni,m)
336  ray = ray + rr(1,i)*ww(ni,l)
337  END DO
338 
339  div(1) = (sum(ff(1,jj(:,m),j) * dw(1,:,l,m)) + &
340  sum(ff(1,jj(:,m),j) * ww(:,l))/ray + &
341  j/ray * sum(ff(4,jj(:,m),j) * ww(:,l)) + &
342  sum(ff(5,jj(:,m),j) * dw(2,:,l,m)))* ray*rj(l,m)
343 
344  IF (j > 0) THEN
345  div(2) = (sum(ff(2,jj(:,m),j) * dw(1,:,l,m)) + &
346  sum(ff(2,jj(:,m),j) * ww(:,l))/ray - &
347  j/ray * sum(ff(3,jj(:,m),j) * ww(:,l)) + &
348  sum(ff(6,jj(:,m),j) * dw(2,:,l,m)))* ray*rj(l,m)
349  ENDIF
350 
351  t = t +(div(1)**2+div(2)**2)
352 
353  ENDDO
354  ENDDO
355 
356  ENDDO
357 
358  t = sqrt(t)
359 
360 
361 END SUBROUTINE nv_1
362 
363 SUBROUTINE nv_0_cn(mesh, ff, p, t)
364  !semi-norme vectoriel H1 pour une representation de Fourier
365  !(v_rc, v_rs, v_zc, v_zs)
366 
367  USE gauss_points
368 
369  IMPLICIT NONE
370 
371  TYPE(mesh_type), TARGET :: mesh
372  REAL(KIND=8), DIMENSION(mesh%np,6), INTENT(IN) :: ff
373  REAL(KIND=8), INTENT(OUT) :: t, p
374 
375  INTEGER :: m, l, i, ni
376  REAL(KIND=8) :: s, rp, rt
377 
378  INTEGER, DIMENSION(:,:), POINTER :: jj
379  INTEGER, POINTER :: me
380  REAL(KIND=8) ,DIMENSION(:,:), POINTER :: rr
381  REAL(KIND=8) :: ray
382 
383  CALL gauss(mesh)
384  jj => mesh%jj
385  me => mesh%me
386  rr=> mesh%rr
387 
388  rp = 0.d0
389  rt = 0.d0
390  s = 0.d0
391  DO m = 1, me
392  DO l = 1, l_g
393  !===Compute radius of Gauss point
394  ray = 0
395  DO ni = 1, n_w; i = jj(ni,m)
396  ray = ray + rr(1,i)*ww(ni,l)
397  END DO
398 
399  rp = rp + sqrt(sum(ff(jj(:,m),1)* ww(:,l))**2 + sum(ff(jj(:,m),5)* ww(:,l))**2)*ray*rj(l,m)
400  rt = rt + abs(sum(ff(jj(:,m),3)* ww(:,l)))*ray*rj(l,m)
401  s = s + ray*rj(l,m)
402 
403  ENDDO
404  ENDDO
405 
406  p = rp/s
407  t = rt/s
408 
409 END SUBROUTINE nv_0_cn
410 
411 END MODULE fem_tn_axi
fem_tn_axi::average
subroutine average(mesh, ff, t)
Definition: fem_tn_axi.f90:138
fem_tn_axi
Definition: fem_tn_axi.f90:4
gauss_points::k_d
integer, public k_d
Definition: associate_gauss.f90:9
gauss_points
Definition: associate_gauss.f90:4
fem_tn_axi::ns_anal_0
subroutine ns_anal_0(mesh, ff, ff_anal, t)
Definition: fem_tn_axi.f90:201
gauss_points::l_g
integer, public l_g
Definition: associate_gauss.f90:9
fem_tn_axi::ns_l1
subroutine ns_l1(mesh, ff, t)
Definition: fem_tn_axi.f90:96
gauss_points::rj
real(kind=8), dimension(:,:), pointer rj
Definition: associate_gauss.f90:14
fem_tn_axi::nv_1
subroutine nv_1(mesh, mod_max, ff, t)
Definition: fem_tn_axi.f90:301
gauss_points::n_w
integer, public n_w
Definition: associate_gauss.f90:9
fem_tn_axi::ns_1
subroutine ns_1(mesh, ff, t)
Definition: fem_tn_axi.f90:249
fem_tn_axi::dot_product
subroutine dot_product(mesh, ff, gg, t)
Definition: fem_tn_axi.f90:9
fem_tn_axi::ns_0
subroutine ns_0(mesh, ff, t)
Definition: fem_tn_axi.f90:52
gauss_points::gauss
subroutine gauss(mesh)
Definition: associate_gauss.f90:24
gauss_points::dw
real(kind=8), dimension(:,:,:,:), pointer dw
Definition: associate_gauss.f90:12
fem_tn_axi::nv_0_cn
subroutine nv_0_cn(mesh, ff, p, t)
Definition: fem_tn_axi.f90:364
gauss_points::ww
real(kind=8), dimension(:,:), pointer ww
Definition: associate_gauss.f90:10