SFEMaNS  version 4.1 (work in progress)
Reference documentation for SFEMaNS
 All Classes Files Functions Variables Groups Pages
maxwell_update_time_with_H.f90
Go to the documentation of this file.
1 !$
2 !Authors Jean-Luc Guermond, Raphael Laguerre, Copyrights 2005
3 !Revised June 2008, Jean-Luc Guermond
4 !Revised Jan/Feb 2009, Caroline Nore, Jean-Luc Guermond, Franky Luddens
5 !
6 MODULE update_maxwell_with_h
7 
8  PUBLIC:: maxwell_decouple_with_h
9  PRIVATE
10  REAL(KIND=8), PARAMETER, PRIVATE :: alpha=0.6d0
11 CONTAINS
12 
13  !------------------------------------------------------------------------------
14  !------------------------------------------------------------------------------
15 
16  SUBROUTINE maxwell_decouple_with_h(comm_one_d, H_mesh, pmag_mesh, phi_mesh, interface_H_phi, &
17  interface_h_mu, hn, bn, phin, hn1, bn1, phin1, vel, stab_in, sigma_in, &
18  r_fourier, index_fourier, mu_h_field, mu_phi, time, dt, rem, list_mode, &
19  h_phi_per, la_h, la_pmag, la_phi, la_mhd, sigma_ns_in, jj_v_to_h)
20  USE def_type_mesh
21  USE chaine_caractere
22  USE solve_petsc
23  USE boundary
24  USE tn_axi
25  USE prep_maill
26  USE dir_nodes_petsc
27  USE st_matrix
28  USE dir_nodes
29  USE my_util
30  USE sft_parallele
31  USE sub_plot
32  USE periodic
33  USE input_data
34  USE verbose
35  IMPLICIT NONE
36  TYPE(mesh_type), INTENT(IN) :: h_mesh, phi_mesh, pmag_mesh
37  TYPE(interface_type), INTENT(IN) :: interface_h_phi, interface_h_mu
38  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
39  REAL(KIND=8), DIMENSION(:,:,:), INTENT(INOUT) :: vel
40  REAL(KIND=8), DIMENSION(H_mesh%np,6,SIZE(list_mode)), INTENT(INOUT) :: hn, hn1
41  REAL(KIND=8), DIMENSION(H_mesh%np,6,SIZE(list_mode)), INTENT(INOUT) :: bn, bn1
42  REAL(KIND=8), DIMENSION(:,:,:), INTENT(INOUT) :: phin, phin1
43  REAL(KIND=8), DIMENSION(3), INTENT(IN) :: stab_in
44  REAL(KIND=8), INTENT(IN) :: r_fourier
45  INTEGER, INTENT(IN) :: index_fourier
46  REAL(KIND=8), INTENT(IN) :: mu_phi, time, dt, rem
47  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: sigma_in, mu_h_field
48  !jan 29/JLG+FL/Forget about it/We replace it by H_p_phi_per/Feb 2 2010
49  TYPE(periodic_type), INTENT(IN) :: h_phi_per
50  !jan 29/JLG+FL/Forget about it/We replace it by H_p_phi_per/Feb 2 2010
51  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: sigma_ns_in
52  INTEGER, DIMENSION(:) , INTENT(IN) :: jj_v_to_h
53  TYPE(petsc_csr_la) :: la_h, la_pmag, la_phi, la_mhd
54  !LC 2016/03/25
55  REAL(KIND=8), DIMENSION(:), ALLOCATABLE, SAVE :: sigma_ns_bar
56  !LC 2016/03/25
57  REAL(KIND=8), DIMENSION(:), ALLOCATABLE, SAVE :: sigma_np
58  REAL(KIND=8), DIMENSION(:), ALLOCATABLE, SAVE :: sigma
59  REAL(KIND=8), SAVE :: sigma_min !FL: not sure we have to save it
60  REAL(KIND=8), DIMENSION(:), POINTER, SAVE :: pmag_bv_d, phi_bv1_d, phi_bv2_d
61  INTEGER, DIMENSION(:), POINTER, SAVE :: pmag_js_d, phi_js_d
62  INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: dirichlet_bdy_h_sides
63  LOGICAL, SAVE :: once=.true.
64  INTEGER, SAVE :: m_max_c
65 !!$ LOGICAL, SAVE :: per = .FALSE.
66  REAL(KIND=8), DIMENSION(3), SAVE :: stab
67  INTEGER, SAVE :: my_petscworld_rank
68  REAL(KIND=8), ALLOCATABLE , DIMENSION(:,:,:), SAVE :: sigma_curl_bdy
69  REAL(KIND=8), DIMENSION(H_mesh%gauss%n_w,H_mesh%me) :: sigma_nj_m
70  REAL(KIND=8), DIMENSION(H_mesh%gauss%l_G*H_mesh%me,6,SIZE(list_mode)) :: sigma_curl
71  REAL(KIND=8), DIMENSION(H_mesh%np,6,SIZE(list_mode)) :: j_over_sigma
72  REAL(KIND=8), DIMENSION(SIZE(Hn,1),6,SIZE(Hn,3)) :: h_ns
73  REAL(KIND=8), DIMENSION(SIZE(Hn,1),2,SIZE(Hn,3)) :: sigma_tot
74  LOGICAL, ALLOCATABLE, DIMENSION(:) :: dir_pmag
75  REAL(KIND=8), DIMENSION(H_mesh%np,6) :: rhs_h
76  REAL(KIND=8), DIMENSION(phi_mesh%np,2) :: rhs_phi
77  REAL(KIND=8), DIMENSION(SIZE(Hn,1),6,SIZE(Hn,3)) :: nl, h_ext
78  REAL(KIND=8), DIMENSION(3) :: temps_par
79  INTEGER, POINTER, DIMENSION(:) :: h_ifrom, pmag_ifrom, phi_ifrom, h_p_phi_ifrom
80  REAL(KIND=8), DIMENSION(phi_mesh%np, 2) :: phin_p1
81  REAL(KIND=8), DIMENSION(H_mesh%np, 6) :: hn_p1
82  INTEGER :: mode, k, i, n, m, ms, code, nj, j
83  INTEGER :: nb_procs, bloc_size, m_max_pad
84  REAL(KIND=8) :: tps, nr_vel, tps_tot, tps_cumul, norm
85  !April 17th, 2008, JLG
86  REAL(KIND=8) :: one_and_half
87  DATA one_and_half/1.5d0/
88  !April 17th, 2008, JLG
89 
90 #include "petsc/finclude/petscsys.h"
91 #include "petsc/finclude/petscmat.h"
92 #include "petsc/finclude/petscksp.h"
93 #include "petsc/finclude/petscvec.h"
94 #include "petsc/finclude/petscvec.h90"
95  petscerrorcode :: ierr
96  mpi_comm, DIMENSION(:), POINTER :: comm_one_d
97  mat, DIMENSION(:), POINTER, SAVE :: h_p_phi_mat1, h_p_phi_mat2
98  mat :: tampon1, tampon2, precond1, precond2
99  ksp, DIMENSION(:), POINTER, SAVE :: h_p_phi_ksp1, h_p_phi_ksp2
100  vec, SAVE :: vx_1, vb_1, vx_1_ghost, vx_2, vb_2, vx_2_ghost
101  !------------------------------END OF DECLARATION--------------------------------------
102 
103  IF (once) THEN
104 
105  once = .false.
106 
107 !!$ IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
108 !!$ per = .TRUE.
109 !!$ ELSE
110 !!$ per = .FALSE.
111 !!$ END IF
112 
113  CALL mpi_comm_rank(petsc_comm_world,my_petscworld_rank,code)
114 
115  CALL create_my_ghost(h_mesh,la_h,h_ifrom)
116  CALL create_my_ghost(pmag_mesh,la_pmag,pmag_ifrom)
117  CALL create_my_ghost(phi_mesh,la_phi,phi_ifrom)
118 
119  n = SIZE(h_ifrom)+SIZE(pmag_ifrom)+SIZE(phi_ifrom)
120  ALLOCATE(h_p_phi_ifrom(n))
121  IF (SIZE(h_ifrom)/=0) THEN
122  h_p_phi_ifrom(1:SIZE(h_ifrom)) = h_ifrom
123  END IF
124  IF (SIZE(pmag_ifrom)/=0) THEN
125  h_p_phi_ifrom(SIZE(h_ifrom)+1:SIZE(h_ifrom)+SIZE(pmag_ifrom)) = pmag_ifrom
126  END IF
127  IF (SIZE(phi_ifrom)/=0) THEN
128  h_p_phi_ifrom(SIZE(h_ifrom)+SIZE(pmag_ifrom)+1:)=phi_ifrom
129  END IF
130 
131  n = 3*h_mesh%dom_np + pmag_mesh%dom_np + phi_mesh%dom_np
132  CALL veccreateghost(comm_one_d(1), n, &
133  petsc_determine, SIZE(h_p_phi_ifrom), h_p_phi_ifrom, vx_1, ierr)
134  CALL vecghostgetlocalform(vx_1, vx_1_ghost, ierr)
135  CALL vecduplicate(vx_1, vb_1, ierr)
136  CALL veccreateghost(comm_one_d(1), n, &
137  petsc_determine, SIZE(h_p_phi_ifrom), h_p_phi_ifrom, vx_2, ierr)
138  CALL vecghostgetlocalform(vx_2, vx_2_ghost, ierr)
139  CALL vecduplicate(vx_2, vb_2, ierr)
140  !------------------------------------------------------------------------------
141 
142  !-------------RESCALING DE SIGMA-----------------------------------------------
143  ALLOCATE(sigma(SIZE(sigma_in)))
144  sigma = sigma_in * rem
145 
146  ! FL, 31/03/11
147  CALL mpi_allreduce(minval(sigma),sigma_min,1,mpi_double_precision, mpi_min,comm_one_d(1), ierr)
148  ! FL, 31/03/11
149 
150  !------------------------------------------------------------------------------
151 
152  !-------------RESCALING DE STAB------------------------------------------------
153  !MARCH, 2010
154  IF (inputs%type_pb=='mhd') THEN
155  ! FL, 31/03/11
156  !stab = stab_in*(1/MINVAL(sigma)+1.d0)
157  stab = stab_in*(1/sigma_min+1.d0)
158  ! FL, 31/03/11
159  ! Velocity assume to be used as reference scale
160 !LC 2016/02/29
161  IF (inputs%if_level_set.AND.inputs%variation_sigma_fluid) THEN
162  stab = stab_in*(1/(minval(inputs%sigma_fluid)*rem)+1.d0)
163  END IF
164 !LC 2016/02/29
165  ELSE
166  nr_vel = norm_sf(comm_one_d, 'L2', h_mesh, list_mode, vel)
167 
168  IF (nr_vel .LE. 1.d-10) THEN
169  ! FL, 31/03/11
170  !stab = stab_in*(1/MINVAL(sigma))
171  stab = stab_in*(1/sigma_min)
172  ! FL, 31/03/11
173  !WRITE(*,*) 'case 1, stab = ',stab
174  ELSE
175  ! FL, 31/03/11
176  !stab = stab_in*(1/MINVAL(sigma)+1.d0)
177  stab = stab_in*(1/sigma_min+1.d0)
178  ! FL, 31/03/11
179  !WRITE(*,*) 'case 2, stab = ',stab
180  ENDIF
181  ! Velocity could be zero in case of Ohmic decay
182  END IF
183  WRITE(*,*) 'stab = ',stab
184  !MARCH, 2010
185 
186  !------------------------------------------------------------------------------
187 
188  !-------------DIMENSIONS-------------------------------------------------------
189  m_max_c = SIZE(list_mode)
190  !------------------------------------------------------------------------------
191 
192  !------------SIGMA IF LEVEL SET------------------------------------------------
193  IF (inputs%if_level_set.AND.inputs%variation_sigma_fluid) THEN
194 !!$ sigma_ns_bar = 0.5d0*MINVAL(inputs%sigma_fluid)* Rem
195 !LC 2016/03/25
196  ALLOCATE(sigma_ns_bar(SIZE(hn,1)))
197 !TEST JLG LC
198 !!$ DO n = 1, SIZE(H_mesh%rr,2)
199 !!$ sigma_ns_bar(n) =0.5d0*MINVAL(inputs%sigma_fluid)*Rem
200 !!$ END DO
201  sigma_ns_bar = sigma_bar_in_fourier_space(h_mesh)*rem
202 !TEST JLG LC
203 !LC 2016/03/25
204  !===check if j=H_mesh%jj(nj,m) is in ns domain or not and define sigma in consequence
205  DO m = 1, h_mesh%me
206  DO nj = 1, h_mesh%gauss%n_w
207  j = h_mesh%jj(nj,m)
208  IF (jj_v_to_h(j) == -1) THEN
209  sigma_nj_m(nj,m) = sigma(m)
210  ELSE
211 !!$ sigma_nj_m(nj,m) = sigma_ns_bar
212 !LC 2016/03/25
213  sigma_nj_m(nj,m) = sigma_ns_bar(j)
214 !LC 2016/03/25
215  END IF
216  END DO
217  END DO
218  ELSE
219  DO m = 1, h_mesh%me
220  sigma_nj_m(:,m) = sigma(m)
221  END DO
222  END IF
223 
224  ALLOCATE(sigma_np(SIZE(hn,1)))
225  sigma_np = 0.d0
226  DO m = 1, h_mesh%me
227  DO nj = 1, h_mesh%gauss%n_w
228  sigma_np(h_mesh%jj(nj,m)) = sigma_nj_m(nj,m)
229  END DO
230  END DO
231  !------------------------------------------------------------------------------
232 
233  !---------------BOUNDARY CONDITIONS FOR pmag----------------------------------
234  ! Creation of Dirichlet boundary conditions for the magnetic pressure
235  ! Only on the boundary that is not periodic...
236  ALLOCATE (dir_pmag(maxval(pmag_mesh%sides)))
237  dir_pmag = .false.
238  DO ms = 1, SIZE(dir_pmag)
239  IF (minval(abs(inputs%list_dirichlet_sides_H-ms)) == 0) THEN
240  dir_pmag(ms) = .true.
241  END IF
242  IF (minval(abs(inputs%list_inter_H_phi-ms)) == 0) THEN
243  dir_pmag(ms) = .true.
244  END IF
245  END DO
246 
247  CALL dirichlet_nodes(pmag_mesh%jjs, pmag_mesh%sides, dir_pmag, pmag_js_d)
248  DEALLOCATE(dir_pmag)
249  ALLOCATE(pmag_bv_d(SIZE(pmag_js_d)))
250  pmag_bv_d = 0.d0
251  ! End creation of Dirichlet boundary conditions for the magnetic pressure
252 
253  !---------------BOUNDARY CONDITIONS FOR Hxn------------------------------------
254  !===Compute sides that are on Dirichlet boundary (H-H_D)xn=0
255  n = 0
256  DO ms = 1, h_mesh%mes
257  IF (minval(abs(h_mesh%sides(ms)-inputs%list_dirichlet_sides_H))/=0) cycle
258  IF (maxval(abs(h_mesh%rr(1,h_mesh%jjs(:,ms)))) .LT.1d-12) cycle
259  n = n + 1
260  END DO
261  ALLOCATE(dirichlet_bdy_h_sides(n))
262  n = 0
263  DO ms = 1, h_mesh%mes
264  IF (minval(abs(h_mesh%sides(ms)-inputs%list_dirichlet_sides_H))/=0) cycle
265  IF (maxval(abs(h_mesh%rr(1,h_mesh%jjs(:,ms)))) .LT.1d-12) cycle
266  n = n + 1
267  dirichlet_bdy_h_sides(n) = ms
268  END DO
269 
270  !---------PREPARE phi_js_D ARRAY FOR POTENTIAL---------------------------------
271  CALL dirichlet_nodes_parallel(phi_mesh, inputs%phi_list_dirichlet_sides, phi_js_d)
272  CALL dirichlet_cavities(comm_one_d(1), interface_h_phi, phi_mesh, phi_js_d)
273  ALLOCATE(phi_bv1_d(SIZE(phi_js_d)), phi_bv2_d(SIZE(phi_js_d)))
274  !------------------------------------------------------------------------------
275 
276  !-------------MATRIX ALLOCATION------------------------------------------------
277  ALLOCATE(h_p_phi_mat1(m_max_c),h_p_phi_ksp1(m_max_c))
278  ALLOCATE(h_p_phi_mat2(m_max_c),h_p_phi_ksp2(m_max_c))
279  IF (SIZE(dirichlet_bdy_h_sides).GE.1) THEN
280  ALLOCATE(sigma_curl_bdy(h_mesh%gauss%l_Gs*SIZE(dirichlet_bdy_h_sides),6,SIZE(list_mode)))
281  ELSE
282  ALLOCATE(sigma_curl_bdy(1,6,SIZE(list_mode)))
283  sigma_curl_bdy = 0.d0
284  END IF
285  !------------------------------------------------------------------------------
286 
287  DO i = 1, m_max_c !Boucle sur les modes
288  mode = list_mode(i)
289 
290  tps = user_time()
291  CALL create_local_petsc_matrix(comm_one_d(1), la_mhd, h_p_phi_mat1(i), clean=.false.)
292  CALL create_local_petsc_matrix(comm_one_d(1), la_mhd, h_p_phi_mat2(i), clean=.false.)
293  IF (i == 1) THEN
294  CALL create_local_petsc_matrix(comm_one_d(1), la_mhd, tampon1, clean=.false.)
295  CALL create_local_petsc_matrix(comm_one_d(1), la_mhd, tampon2, clean=.false.)
296  CALL create_local_petsc_matrix(comm_one_d(1), la_mhd, precond1, clean=.false.)
297  CALL create_local_petsc_matrix(comm_one_d(1), la_mhd, precond2, clean=.false.)
298  END IF
299 
300  tps = user_time() - tps
301 !!$ WRITE(*,*) ' Tps create_local_petsc_matrix', tps
302 
303  tps = user_time()
304  CALL mat_h_p_phi_maxwell(h_mesh,pmag_mesh,phi_mesh,interface_h_phi, &
305  mode,mu_h_field, mu_phi, one_and_half/dt, stab, r_fourier, index_fourier, &
306  la_h, la_pmag, la_phi, h_p_phi_mat1(i), h_p_phi_mat2(i), sigma_np)
307 
308  tps = user_time() - tps
309 !!$ WRITE(*,*) ' Tps mat_H_p_phi_maxwell', tps
310 
311  !Take care of discontinuous mu
312  tps = user_time()
313  CALL mat_maxwell_mu(h_mesh, interface_h_mu, mode, stab, &
314  mu_h_field, sigma, la_h, h_p_phi_mat1(i), h_p_phi_mat2(i))
315  tps = user_time() - tps
316 !!$ WRITE(*,*) ' Tps mat_maxwell_mu', tps
317  !JLG, FL, Feb 10, 2011
318 
319  tps = user_time()
320  CALL mat_dirichlet_maxwell(h_mesh, dirichlet_bdy_h_sides, &
321  mode, stab, la_h, h_p_phi_mat1(i), h_p_phi_mat2(i), sigma_np)
322 
323 !!$ tps = user_time() - tps
324 !!$ WRITE(*,*) ' Tps mat_dirichlet_maxwell', tps
325 
326 !!$ IF (per) THEN
327  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
328  CALL periodic_matrix_petsc(h_phi_per%n_bord, h_phi_per%list, &
329  h_phi_per%perlist, h_p_phi_mat1(i), la_mhd)
330  CALL periodic_matrix_petsc(h_phi_per%n_bord, h_phi_per%list, &
331  h_phi_per%perlist, h_p_phi_mat2(i), la_mhd)
332  END IF
333  tps = user_time()
334  CALL dirichlet_m_parallel(h_p_phi_mat1(i),la_pmag%loc_to_glob(1,pmag_js_d))
335  CALL dirichlet_m_parallel(h_p_phi_mat1(i),la_phi%loc_to_glob(1,phi_js_d))
336  CALL dirichlet_m_parallel(h_p_phi_mat2(i),la_pmag%loc_to_glob(1,pmag_js_d))
337  CALL dirichlet_m_parallel(h_p_phi_mat2(i),la_phi%loc_to_glob(1,phi_js_d))
338  tps = user_time() - tps
339 
340 !!$ WRITE(*,*) ' Tps Dirichlet_M_parallel', tps
341 
342  tps = user_time()
343  CALL init_solver(inputs%my_par_H_p_phi,h_p_phi_ksp1(i),h_p_phi_mat1(i),comm_one_d(1),&
344  solver=inputs%my_par_H_p_phi%solver,precond=inputs%my_par_H_p_phi%precond)
345  CALL init_solver(inputs%my_par_H_p_phi,h_p_phi_ksp2(i),h_p_phi_mat2(i),comm_one_d(1),&
346  solver=inputs%my_par_H_p_phi%solver,precond=inputs%my_par_H_p_phi%precond)
347  tps = user_time() - tps
348 !!$ WRITE(*,*) ' Tps init_solver', tps
349 
350 !!$ !==================TEST===================
351  CALL matdestroy(h_p_phi_mat1(i),ierr)
352  CALL matdestroy(h_p_phi_mat2(i),ierr)
353 
354  ENDDO
355 
356  !------------------------------------------------------------------------------
357  ENDIF ! fin du once
358  tps_tot = user_time()
359  tps_cumul = 0
360  CALL mpi_comm_rank(petsc_comm_world, my_petscworld_rank, code)
361 
362  !-------------TRANSPORT TERM---------------------------------------------------
363  tps = user_time()
364  nr_vel = norm_sf(comm_one_d, 'L2', h_mesh, list_mode, vel)
365  h_ext = 2*hn - hn1
366 
367  IF (nr_vel .LE. 1.d-10) THEN
368  nl = 0.d0
369  ELSE
370  CALL mpi_comm_size(comm_one_d(2), nb_procs, code)
371  bloc_size = SIZE(vel,1)/nb_procs+1
372  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
373  CALL fft_par_cross_prod_dcl(comm_one_d(2), vel, h_ext, nl, nb_procs, bloc_size, m_max_pad, temps_par)
374  ENDIF
375 
376  IF (inputs%if_level_set.AND.inputs%variation_sigma_fluid) THEN
377  h_ns = 0.d0
378  sigma_tot = 0.d0
379  DO m = 1, h_mesh%me
380  DO nj = 1, h_mesh%gauss%n_w
381  j = h_mesh%jj(nj,m)
382  !Check if node is in Navier-Stokes domain(s)
383  IF (jj_v_to_h(j) /= -1) THEN
384  h_ns(j,:,:) = 2*hn(j,:,:)- hn1(j,:,:)
385  sigma_tot(j,:,:) = sigma_ns_in(jj_v_to_h(j),:,:) * rem
386  ELSE
387  DO i = 1, SIZE(list_mode)
388  mode = list_mode(i)
389  IF (mode==0) THEN
390  sigma_tot(j,1,i) = sigma(m)
391  END IF
392  END DO
393  END IF
394  END DO
395  END DO
396 
397  !===Compute (1/sigma_ns_bar - 1/sigma)*CURL(H_ns) in fluid domain and 0 elsewhere
398  CALL smb_sigma_prod_curl(comm_one_d(2), h_mesh, list_mode, h_ns, &
399  sigma_ns_bar, sigma_tot, sigma_curl)
400  IF (SIZE(dirichlet_bdy_h_sides).GE.1) THEN
401  CALL smb_sigma_prod_curl_bdy(comm_one_d(2), h_mesh, dirichlet_bdy_h_sides, list_mode, h_ns, &
402  sigma_ns_bar, sigma_tot, sigma_curl_bdy)
403  ELSE
404  sigma_curl_bdy = 0.d0
405  END IF
406  !===Compute J/sigma
407  CALL smb_current_over_sigma(comm_one_d(2), h_mesh, list_mode, &
408  mu_h_field, mu_phi, sigma_tot, time, j_over_sigma)
409  ELSE
410  sigma_curl = 0.d0
411  sigma_curl_bdy = 0.d0
412  j_over_sigma = 0.d0
413  END IF
414 
415  tps = user_time() - tps; tps_cumul=tps_cumul+tps
416  !WRITE(*,*) ' Tps NLS_SFT Maxwell', tps
417  !------------------------------------------------------------------------------
418 
419  !-------------SOLUTION OF LINEAR SYSTEMS---------------------------------------
420  DO i = 1, m_max_c
421 
422  mode = list_mode(i)
423 
424  !-------------SOURCES TERMS----------------------------------------------------
425  tps = user_time()
426  DO k = 1, 6
427  !rhs_H (:,k) = mu_H_field*(4*Hn(:,k,i)-Hn1(:,k,i))/(2*dt)
428  rhs_h(:,k) = (4*bn(:,k,i)-bn1(:,k,i))/(2*dt)
429  END DO
430  DO k = 1, 2
431  rhs_phi(:,k) = mu_phi*(4*phin(:,k,i)-phin1(:,k,i))/(2*dt)
432  END DO
433  !-------------Integration by parts of the scalar potential------------------
434  CALL courant_int_by_parts(h_mesh,phi_mesh,interface_h_phi,sigma,mu_phi,mu_h_field,time,mode, &
435  rhs_h, nl(:,:,i), la_h, la_phi, vb_1, vb_2, h_ext(:,:,i), &
436  sigma_curl(:,:,i), j_over_sigma(:,:,i))
437 !!$ ! Feb 2010, JLG + FL
438 !!$ CALL courant_int_by_parts(H_mesh,phi_mesh,interface_H_phi,sigma,mu_phi,mu_H_field,time,mode, &
439 !!$ rhs_H, NL(:,:,i), LA_H, LA_phi, vb_1, vb_2, H_ext(:,:,i))
440 
441  !-------------Integration by parts of the scalar potential------------------
442 
443  tps = user_time() - tps; tps_cumul=tps_cumul+tps
444  !WRITE(*,*) ' Tps courant', tps
445  !Takes care of discontinuous mu
446  tps = user_time()
447  CALL courant_mu(h_mesh, interface_h_mu,sigma, mu_h_field, time,mode, nl(:,:,i), &
448  la_h, vb_1, vb_2, h_ext(:,:,i))
449  tps = user_time() - tps; tps_cumul=tps_cumul+tps
450 !!$ WRITE(*,*) ' Tps courant_mu', tps
451 
452  !JLG, FL, Feb 10, 2011
453  !Take care of Dirichlet conditions on H (H x n = Hd x n)
454  !CALL rhs_dirichlet(H_mesh, Dirichlet_bdy_H_sides, &
455  ! sigma, mu_H_field, time, mode, NL, i, stab, LA_H, vb_1,vb_2,H_ext)
456  CALL rhs_dirichlet(h_mesh, dirichlet_bdy_h_sides, &
457  sigma, mu_h_field, time, mode, nl(:,:,i), stab, la_h, vb_1,vb_2, h_ext(:,:,i), &
458  j_over_sigma(:,:,i), sigma_curl_bdy(:,:,i))
459  !------------------------------------------------------------------------------
460 
461  !-------------INTERFACE INTEGRAL-----------------------------------------------
462  tps = user_time()
463  CALL surf_int(h_mesh,phi_mesh,interface_h_phi,interface_h_mu,inputs%list_dirichlet_sides_H, &
464  sigma,mu_phi,mu_h_field, time,mode,la_h, la_phi,vb_1,vb_2, r_fourier, index_fourier)
465  tps = user_time() - tps; tps_cumul=tps_cumul+tps
466 !!$ WRITE(*,*) ' Tps surf_int', tps
467  !------------------------------------------------------------------------------
468 
469  !---------------------PERIODIC-------------------
470 !!$ IF (per) THEN
471  IF (inputs%my_periodic%nb_periodic_pairs/=0) THEN
472  CALL periodic_rhs_petsc(h_phi_per%n_bord, h_phi_per%list, h_phi_per%perlist, vb_1, la_mhd)
473  CALL periodic_rhs_petsc(h_phi_per%n_bord, h_phi_per%list, h_phi_per%perlist, vb_2, la_mhd)
474  END IF
475 
476  !-------------DIRICHLET BOUNDARY CONDITIONS-------------------------------------
477  tps = user_time()
478  CALL dirichlet_rhs(la_pmag%loc_to_glob(1,pmag_js_d)-1, pmag_bv_d,vb_1)
479  CALL dirichlet_rhs(la_pmag%loc_to_glob(1,pmag_js_d)-1, pmag_bv_d,vb_2)
480  IF (SIZE(phi_js_d)>0) THEN
481  phi_bv1_d = phiexact(1,phi_mesh%rr(1:2,phi_js_d), mode, mu_phi, time)
482  phi_bv2_d = phiexact(2,phi_mesh%rr(:,phi_js_d), mode, mu_phi, time)
483  ELSE
484  phi_bv1_d = 0.d0
485  phi_bv2_d = 0.d0
486  END IF
487 
488  CALL dirichlet_rhs(la_phi%loc_to_glob(1,phi_js_d)-1, phi_bv1_d, vb_1)
489  CALL dirichlet_rhs(la_phi%loc_to_glob(1,phi_js_d)-1, phi_bv2_d, vb_2)
490  !CALL dirichlet_rhs(LA_phi%loc_to_glob(1,phi_js_D)-1,&
491  ! Phiexact(1,phi_mesh%rr(:,phi_js_D), mode, mu_phi, time),vb_1)
492  !CALL dirichlet_rhs(LA_phi%loc_to_glob(1,phi_js_D)-1,&
493  ! Phiexact(2,phi_mesh%rr(:,phi_js_D), mode, mu_phi, time),vb_2)
494 
495  tps = user_time() - tps; tps_cumul=tps_cumul+tps
496 !!$ WRITE(*,*) ' Tps bcs', tps
497  !-------------------------------------------------------------------------------
498 
499  !-------------SOLVING THE LINEAR SYSTEMS----------------------------------------
500  IF (inputs%my_par_H_p_phi%verbose .AND. (i==1)) WRITE(*,*) 'start solving'
501  tps = user_time()
502 
503  CALL solver(h_p_phi_ksp1(i),vb_1,vx_1,reinit=.false.,verbose=inputs%my_par_H_p_phi%verbose)
504 
505  CALL vecghostupdatebegin(vx_1,insert_values,scatter_forward,ierr)
506  CALL vecghostupdateend(vx_1,insert_values,scatter_forward,ierr)
507  IF (h_mesh%me/=0) THEN
508  CALL extract(vx_1_ghost,1,1,la_mhd,hn_p1(:,1))
509  CALL extract(vx_1_ghost,2,2,la_mhd,hn_p1(:,4))
510  CALL extract(vx_1_ghost,3,3,la_mhd,hn_p1(:,5))
511  END IF
512  IF (phi_mesh%me/=0) THEN
513  CALL extract(vx_1_ghost,5,5,la_mhd,phin_p1(:,1))
514  END IF
515 
516  CALL solver(h_p_phi_ksp2(i),vb_2,vx_2,reinit=.false.,verbose=inputs%my_par_H_p_phi%verbose)
517  CALL vecghostupdatebegin(vx_2,insert_values,scatter_forward,ierr)
518  CALL vecghostupdateend(vx_2,insert_values,scatter_forward,ierr)
519  IF (h_mesh%me/=0) THEN
520  CALL extract(vx_2_ghost,1,1,la_mhd,hn_p1(:,2))
521  CALL extract(vx_2_ghost,2,2,la_mhd,hn_p1(:,3))
522  CALL extract(vx_2_ghost,3,3,la_mhd,hn_p1(:,6))
523  END IF
524  IF (phi_mesh%me/=0) THEN
525  CALL extract(vx_2_ghost,5,5,la_mhd,phin_p1(:,2))
526  END IF
527 
528  tps = user_time() - tps; tps_cumul=tps_cumul+tps
529  !WRITE(*,*) ' Tps solve Maxwell', tps
530  !-------------------------------------------------------------------------------
531 
532 
533  !-------------UPDATE------------------------------------------------------------
534  !JLG AR, Dec 18 2008
535  IF (mode==0) THEN
536  IF (h_mesh%me /=0) THEN
537  hn_p1(:,2) = 0.d0
538  hn_p1(:,4) = 0.d0
539  hn_p1(:,6) = 0.d0
540  END IF
541  IF (phi_mesh%me /=0 ) THEN
542  phin_p1(:,2) = 0.d0
543  END IF
544  END IF
545  !JLG AR, Dec 18 2008
546 
547  tps = user_time()
548  IF (h_mesh%me /=0) THEN
549  hn1(:,:,i) = hn(:,:,i)
550 
551  hn(:,1,i) = hn_p1(:,1)
552  hn(:,4,i) = hn_p1(:,4)
553  hn(:,5,i) = hn_p1(:,5)
554 
555  hn(:,2,i) = hn_p1(:,2)
556  hn(:,3,i) = hn_p1(:,3)
557  hn(:,6,i) = hn_p1(:,6)
558 
559  DO k = 1, 6
560  bn1(:,k,i) = bn(:,k,i)
561  bn(:,k,i) = mu_h_field*hn(:,k,i)
562  END DO
563 
564  END IF
565 
566  IF (phi_mesh%me /= 0) THEN
567  phin1(:,:,i) = phin(:,:,i)
568 
569  phin(:,1,i) = phin_p1(:,1)
570 
571  phin(:,2,i) = phin_p1(:,2)
572  END IF
573  tps = user_time() - tps; tps_cumul=tps_cumul+tps
574  !WRITE(*,*) ' Tps update', tps
575  !------------------------------------------------------------------------------
576 
577  ENDDO
578 
579  !===Verbose divergence of velocity
580  IF (inputs%verbose_divergence) THEN
581  norm = norm_sf(comm_one_d, 'L2', h_mesh, list_mode, bn)
582  talk_to_me%div_B_L2 = norm_sf(comm_one_d, 'div', h_mesh, list_mode, bn)/norm
583  talk_to_me%time=time
584  END IF
585 
586  tps_tot = user_time() - tps_tot
587 !!$ WRITE(*,'(A,2(f13.3,2x),10(I3,x))') ' Tps boucle en temps Maxwell', tps_tot, tps_cumul, list_mode
588 !!$ WRITE(*,*) ' TIME = ', time, '========================================'
589 
590  END SUBROUTINE maxwell_decouple_with_h
591 
592  SUBROUTINE mat_h_p_phi_maxwell(H_mesh, pmag_mesh, phi_mesh, interface_H_phi, &
593  mode, mu_h_field, mu_phi, c_mass, stab, r_fourier, index_fourier, &
594  la_h, la_pmag, la_phi, h_p_phi_mat1, h_p_phi_mat2, sigma_np)
595  USE def_type_mesh
596  USE dir_nodes
597  USE gauss_points
598  USE boundary
599  IMPLICIT NONE
600  TYPE(mesh_type), INTENT(IN) :: h_mesh
601  TYPE(mesh_type), INTENT(IN) :: pmag_mesh
602  TYPE(mesh_type), INTENT(IN) :: phi_mesh
603  TYPE(interface_type), INTENT(IN) :: interface_h_phi
604  INTEGER, INTENT(IN) :: mode
605  REAL(KIND=8), INTENT(IN) :: mu_phi, c_mass
606  REAL(KIND=8), DIMENSION(3), INTENT(IN) :: stab
607  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: mu_h_field
608  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: sigma_np
609  REAL(KIND=8), OPTIONAL :: r_fourier
610  INTEGER, OPTIONAL :: index_fourier
611  REAL(KIND=8), DIMENSION(phi_mesh%gauss%n_ws,phi_mesh%gauss%l_Gs) :: w_cs
612  REAL(KIND=8), DIMENSION(2, H_mesh%gauss%n_w, phi_mesh%gauss%l_Gs, H_mesh%mes) :: dw_cs
613  INTEGER :: m, l, ms, ls, ni, nj, k, i, j, &
614  n_ws1, n_ws2, n_w2, n_w1, m1, m2, ki, kj,ib,jb, ms1, ms2
615  REAL(KIND=8) :: x, y, hm1, stab_div, stab_colle_h_phi
616  REAL(KIND=8) :: ray, error
617  LOGICAL :: mark=.false.
618  REAL(KIND=8), DIMENSION(3,H_mesh%gauss%n_w,pmag_mesh%gauss%n_w) :: thpmag
619  REAL(KIND=8), DIMENSION(pmag_mesh%gauss%n_w,pmag_mesh%gauss%n_w) :: tpmag
620  REAL(KIND=8), DIMENSION(9,H_mesh%gauss%n_w,H_mesh%gauss%n_w) :: th
621  REAL(KIND=8), DIMENSION(phi_mesh%gauss%n_w,phi_mesh%gauss%n_w):: tphi
622 
623  !MATRICES POUR LES TERMES DE VOLUMES c_mass*mu_H*H + Rot((1/sigma)Rot(H)) - Grad(Div(H))
624  ! -c_mass*mu_phi*Lap(Phi)
625  !========================================================================
626  !Le probleme est decouple en deux sous groupes de variables :
627  !H1, H4, H5 et Phi1 d'une part et H2, H3, H6 et Phi2 d'autre part.
628  !Les matrices (symetriques sans terme de bord) s'ecrivent :
629 
630  !MATRICE 1 ::
631  ! (------------------------------)
632  ! ( TH1 | TH2 | TH3 | | ) H1
633  ! ( | TH4 | TH5 | | ) H4
634  ! ( | TH6 | | ) H5
635  ! ( | Tpmag | ) P1
636  ! ( |TPhi) Phi1
637  ! (------------------------------)
638 
639  !MATRICE 2 (TH2 => TH8 et TH5 => TH9::
640  ! (------------------------)
641  ! ( TH1 | TH8 | TH3 | ) H2
642  ! ( | TH4 | TH9 | ) H3
643  ! ( | TH6 | ) H6
644  ! ( | TPhi ) Phi2
645  ! (------------------------)
646  !=========================================================================
647 
648 
649  REAL(KIND=8), DIMENSION(9,H_mesh%gauss%n_w,H_mesh%gauss%n_w) :: hsij
650  REAL(KIND=8), DIMENSION(phi_mesh%gauss%n_w,phi_mesh%gauss%n_w) :: phisij
651  REAL(KIND=8), DIMENSION(6,phi_mesh%gauss%n_w,phi_mesh%gauss%n_w) :: sij
652 
653  ! MATRICES POUR LES TERMES DE BORDS Hsij et Phisij
654  !=================================================
655  ! (--------------------------------------------------------------------)
656  ! ( Hsij(1) | Hsij(2) | Hsij(4) || Sij(1) )
657  ! ( Hsij(1) | Hsij(3) | Hsij(4) || Sij(2) )
658  ! (--------------------------------------------------------------------)
659  ! ( | Hsij(5) | || Sij(3) )
660  ! ( | Hsij(5) | || Sij(4) )
661  ! (--------------------------------------------------------------------)
662  ! ( Hsij(7) | Hsij(9) | Hsij(6) || Sij(5) )
663  ! ( Hsij(7) | Hsij(8) | Hsij(6) || Sij(6) )
664  ! (====================================================================)
665  ! ( Sij'(1) | Sij'(3) | Sij'(5) || Phisij )
666  ! ( Sij'(2) | Sij'(4) | Sij'(6) || Phisij )
667  ! (------------------------------------------------------------------- )
668  !
669  ! L'autre partie des termes croises est la symetrique de la premiere
670  ! juste apres le calcul du terme de bord dissymetrique
671 
672  !fonctions de forme propres a H_mesh
673  REAL(KIND=8), DIMENSION(:,:), POINTER :: ww_h
674  !derivees des fonctions de forme propres a H_mesh
675  REAL(KIND=8), DIMENSION(:,:,:,:), POINTER :: dw_h
676  !jacobien pour H
677  REAL(KIND=8), DIMENSION(:,:), POINTER :: rj_h
678  !fonctions de forme propres a phi_mesh
679  REAL(KIND=8), DIMENSION(:,:), POINTER :: ww_phi
680  !derivees des fonctions de forme propres a phi_mesh
681  REAL(KIND=8), DIMENSION(:,:,:,:), POINTER :: dw_phi
682  REAL(KIND=8), DIMENSION(2,H_mesh%gauss%n_w,H_mesh%gauss%l_G) :: dwp
683  REAL(KIND=8), DIMENSION(H_mesh%gauss%n_w,H_mesh%gauss%l_G) :: wwp
684  !jacobien pour phi
685  REAL(KIND=8), DIMENSION(:,:), POINTER :: rj_phi
686 
687  REAL(KIND=8), DIMENSION(2,phi_mesh%gauss%l_Gs) :: gauss1, gauss2
688  INTEGER :: ls1, ls2
689  REAL(KIND=8) :: ref, diff, mu_h, c_mu_h, c_mu_phi, muhl, &
690  dzmuhl, drmuhl, c_div, hloc, viscolm, xij, eps
691  !June 8 2008
692  REAL(KIND=8) :: c_sym=.0d0 ! Symmetrization of the bilinear form
693  !June 8 2008
694  !June 2009, JLG, CN, Normalization
695  REAL(KIND=8) :: c_lap
696  !June 2009, JLG, CN
697 !!$ FL + CN 22/03/2013
698 !!$ REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: mat_loc1, mat_loc2
699 !!$ INTEGER , DIMENSION(:), ALLOCATABLE :: idxn, jdxn
700  REAL(KIND=8), DIMENSION(3*H_mesh%gauss%n_w+pmag_mesh%gauss%n_w+ & phi_mesh%gauss%n_w , 3*H_mesh%gauss%n_w+pmag_mesh%gauss%n_w+ & phi_mesh%gauss%n_w) :: mat_loc1, mat_loc2
701  INTEGER , DIMENSION(3*H_mesh%gauss%n_w+pmag_mesh%gauss%n_w+ & phi_mesh%gauss%n_w) :: idxn, jdxn
702 !!$ FL + CN 22/03/2013
703  TYPE(petsc_csr_la) :: la_h, la_pmag, la_phi
704  INTEGER :: n_wpmag, n_wh, n_wphi, ix, jx
705 !LC 2016/03/25
706  REAL(KIND=8) :: sigma_np_gauss
707 !LC 2016/03/25
708 
709 #include "petsc/finclude/petsc.h"
710  mat :: h_p_phi_mat1, h_p_phi_mat2
711  petscerrorcode :: ierr
712 
713  CALL matzeroentries(h_p_phi_mat1, ierr)
714  CALL matzeroentries(h_p_phi_mat2, ierr)
715  CALL matsetoption(h_p_phi_mat1, mat_row_oriented, petsc_false, ierr)
716  CALL matsetoption(h_p_phi_mat2, mat_row_oriented, petsc_false, ierr)
717 
718  !June 2009, JLG, CN, Normalization
719  c_lap = .1d0
720  stab_colle_h_phi = stab(2)
721  stab_div = stab(1)
722  !Jan 2010, JLG, CN, Normalization,
723 
724  c_mu_phi = c_mass*mu_phi
725 
726  ww_h => h_mesh%gauss%ww
727  dw_h => h_mesh%gauss%dw
728  rj_h => h_mesh%gauss%rj
729  ww_phi => phi_mesh%gauss%ww
730  dw_phi => phi_mesh%gauss%dw
731  rj_phi => phi_mesh%gauss%rj
732 
733  n_wh = h_mesh%gauss%n_w
734  n_wpmag = pmag_mesh%gauss%n_w
735  n_wphi = phi_mesh%gauss%n_w
736 
737  !==Block on H
738 !!$ ALLOCATE(mat_loc1(3*n_wH+n_wpmag+n_wphi,3*n_wH+n_wpmag+n_wphi))
739 !!$ ALLOCATE(mat_loc2(3*n_wH+n_wpmag+n_wphi,3*n_wH+n_wpmag+n_wphi))
740 !!$ ALLOCATE(jdxn(3*n_wH+n_wpmag+n_wphi),idxn(3*n_wH+n_wpmag+n_wphi))
741  DO m = 1, h_mesh%me
742 
743  th = 0.d0
744 
745  DO l = 1, h_mesh%gauss%l_G
746  hloc = sqrt(sum(h_mesh%gauss%rj(:,m)))**(2*alpha)
747  !===Compute radius of Gauss point
748  !Feb 8 2007, muhl
749  muhl = sum(mu_h_field(h_mesh%jj(:,m))*ww_h(:,l))
750  drmuhl = sum(mu_h_field(h_mesh%jj(:,m))*dw_h(1,:,l,m))
751  dzmuhl = sum(mu_h_field(h_mesh%jj(:,m))*dw_h(2,:,l,m))
752  c_mu_h = c_mass*muhl
753  !Feb 8 2007, muhl
754 !LC 2016/03/28
755  sigma_np_gauss = sum(sigma_np(h_mesh%jj(:,m))*ww_h(:,l))
756 !LC 2016/03/28
757  !June 7 2008, Normalization, JLG, FL, May, 28, 2009
758  c_div = stab_div*hloc
759  !c_div = stab_div*hloc/muhl
760  !c_div = stab_div*hloc/muhl**2
761  !c_div = stab_div
762  !c_div = stab_div/muhl
763  !c_div = stab_div/muhl**2
764  !June 7 2008, Normalization
765 
766  ray = 0
767  DO ni = 1, h_mesh%gauss%n_w; i = h_mesh%jj(ni,m)
768  ray = ray + h_mesh%rr(1,i)*ww_h(ni,l)
769  END DO
770 
771  DO ni = 1, h_mesh%gauss%n_w
772  DO nj = 1, h_mesh%gauss%n_w
773  j = h_mesh%jj(nj,m)
774 
775 ! TEST
776  ! mu_H * <bi,bj> + <Div bi,Div bj> + <(1/sigma) Rot bi,Rot bj>
777  th(1,ni,nj) = th(1,ni,nj) + rj_h(l,m) * ray* ( &
778 !DCQ + JLG (Nov 13 2013). Mass integration done same way on LHS and RHS.
779 ! c_mu_H*ww_H(ni,l)*ww_H(nj,l) &
780  c_mass*mu_h_field(j)*ww_h(nj,l)*ww_h(ni,l) &
781  + (dw_h(2,ni,l,m)*dw_h(2,nj,l,m) + mode**2/ray**2*ww_h(ni,l)*ww_h(nj,l))/sigma_np_gauss &
782  !DIVERGENCE, June 8 2008
783  + c_div*(muhl*(ww_h(ni,l)/ray+dw_h(1,ni,l,m)) + ww_h(ni,l)*drmuhl) &
784  *(muhl*(ww_h(nj,l)/ray+dw_h(1,nj,l,m)) + ww_h(nj,l)*drmuhl))
785  !+ stab_div*(ww_H(ni,l)*ww_H(nj,l)/ray**2+dw_H(1,ni,l,m)*dw_H(1,nj,l,m) &
786  !+ 1/ray*(ww_H(ni,l)*dw_H(1,nj,l,m)+ww_H(nj,l)*dw_H(1,ni,l,m))))
787  !
788 
789  th(2,ni,nj) = th(2,ni,nj)+ rj_h(l,m) * ray* ( &
790  mode/ray**2 * ww_h(ni,l)*(ww_h(nj,l)+ray*dw_h(1,nj,l,m))/sigma_np_gauss &
791  !DIVERGENCE , June 8 2008
792  + c_div*mode/ray*(muhl*(ww_h(ni,l)/ray+dw_h(1,ni,l,m)) &
793  + ww_h(ni,l)*drmuhl)*muhl*ww_h(nj,l))
794  !+ stab_div*mode/ray*(ww_H(ni,l)/ray+dw_H(1,ni,l,m))*ww_H(nj,l))
795  !
796 
797  th(8,ni,nj) = th(8,ni,nj)+ rj_h(l,m) * ray* ( &
798  - mode/ray**2 * ww_h(ni,l)*(ww_h(nj,l)+ray*dw_h(1,nj,l,m))/sigma_np_gauss &
799  !DIVERGENCE, June 8 2008
800  - c_div*mode/ray*(muhl*(ww_h(ni,l)/ray+dw_h(1,ni,l,m)) &
801  + ww_h(ni,l)*drmuhl)*muhl*ww_h(nj,l))
802  !-stab_div*mode/ray*(ww_H(ni,l)/ray+dw_H(1,ni,l,m))*ww_H(nj,l))
803  !
804 
805  th(3,ni,nj) = th(3,ni,nj)+ rj_h(l,m) * ray* ( &
806  - dw_h(2,ni,l,m)*dw_h(1,nj,l,m)/sigma_np_gauss &
807  !DIVERGENCE, June 8 2008
808  + c_div*(muhl*(ww_h(ni,l)/ray+dw_h(1,ni,l,m)) + ww_h(ni,l)*drmuhl)*&
809  (muhl*dw_h(2,nj,l,m) + ww_h(nj,l)*dzmuhl))
810  !+ stab_div*(ww_H(ni,l)/ray+dw_H(1,ni,l,m))*dw_H(2,nj,l,m))
811  !
812 
813  th(4,ni,nj) = th(4,ni,nj) + rj_h(l,m) * ray* ( &
814 ! c_mu_H*ww_H(ni,l)*ww_H(nj,l) &
815 !DCQ + JLG (Nov 13 2013). Mass integration done same way on LHS and RHS.
816  c_mass*mu_h_field(j)*ww_h(nj,l)*ww_h(ni,l) &
817  + (dw_h(2,ni,l,m)*dw_h(2,nj,l,m) &
818  + 1/ray**2 *(ww_h(ni,l)+ray*dw_h(1,ni,l,m))*(ww_h(nj,l)&
819  +ray*dw_h(1,nj,l,m)))/sigma_np_gauss &
820  !DIVERGENCE, June 8 2008
821  +c_div*muhl**2*mode**2/ray**2*ww_h(ni,l)*ww_h(nj,l))
822  !+stab_div*mode**2/ray**2*ww_H(ni,l)*ww_H(nj,l))
823  !
824 
825  th(5,ni,nj) = th(5,ni,nj) + rj_h(l,m) * ray* (&
826  + mode/ray*dw_h(2,ni,l,m)*ww_h(nj,l)/sigma_np_gauss &
827  !DIVERGENCE, June 8 2008
828  +c_div*mode/ray*muhl*ww_h(ni,l)*(muhl*dw_h(2,nj,l,m) + ww_h(nj,l)*dzmuhl))
829  !+stab_div*mode/ray*ww_H(ni,l)*dw_H(2,nj,l,m))
830  !
831 
832  th(9,ni,nj) = th(9,ni,nj) + rj_h(l,m) * ray* (&
833  - mode/ray*dw_h(2,ni,l,m)*ww_h(nj,l)/sigma_np_gauss &
834  !DIVERGENCE, June 8 2008
835  - c_div*mode/ray*muhl*ww_h(ni,l)*(muhl*dw_h(2,nj,l,m) + ww_h(nj,l)*dzmuhl))
836  !- stab_div*mode/ray*ww_H(ni,l)*dw_H(2,nj,l,m))
837  !
838 
839  th(6,ni,nj) = th(6,ni,nj) + rj_h(l,m) * ray* ( &
840 ! c_mu_H*ww_H(ni,l)*ww_H(nj,l) &
841 !DCQ + JLG (Nov 13 2013). Mass integration done same way on LHS and RHS.
842  c_mass*mu_h_field(j)*ww_h(nj,l)*ww_h(ni,l) &
843  + (mode**2/ray**2*ww_h(ni,l)*ww_h(nj,l) + dw_h(1,ni,l,m)*dw_h(1,nj,l,m))/sigma_np_gauss &
844  !DIVERGENCE, June 8 2008
845  + c_div*(muhl*dw_h(2,ni,l,m) + ww_h(ni,l)*dzmuhl) &
846  *(muhl*dw_h(2,nj,l,m) + ww_h(nj,l)*dzmuhl))
847  !+ stab_div*dw_H(2,ni,l,m)*dw_H(2,nj,l,m))
848  !
849 ! TEST
850 
851 !!$ ! mu_H * <bi,bj> + <Div bi,Div bj> + <(1/sigma) Rot bi,Rot bj>
852 !!$ TH(1,ni,nj) = TH(1,ni,nj) + rj_H(l,m) * ray* ( &
853 !!$ c_mu_H*ww_H(ni,l)*ww_H(nj,l) &
854 !!$ + (dw_H(2,ni,l,m)*dw_H(2,nj,l,m) + mode**2/ray**2*ww_H(ni,l)*ww_H(nj,l))/sigma(m) &
855 !!$ !DIVERGENCE, June 8 2008
856 !!$ + c_div*(muhl*(ww_H(ni,l)/ray+dw_H(1,ni,l,m)) + ww_H(ni,l)*drmuhl) &
857 !!$ *(muhl*(ww_H(nj,l)/ray+dw_H(1,nj,l,m)) + ww_H(nj,l)*drmuhl))
858 !!$ !+ stab_div*(ww_H(ni,l)*ww_H(nj,l)/ray**2+dw_H(1,ni,l,m)*dw_H(1,nj,l,m) &
859 !!$ !+ 1/ray*(ww_H(ni,l)*dw_H(1,nj,l,m)+ww_H(nj,l)*dw_H(1,ni,l,m))))
860 !!$ !
861 !!$
862 !!$ TH(2,ni,nj) = TH(2,ni,nj)+ rj_H(l,m) * ray* ( &
863 !!$ mode/ray**2 * ww_H(ni,l)*(ww_H(nj,l)+ray*dw_H(1,nj,l,m))/sigma(m) &
864 !!$ !DIVERGENCE , June 8 2008
865 !!$ + c_div*mode/ray*(muhl*(ww_H(ni,l)/ray+dw_H(1,ni,l,m)) &
866 !!$ + ww_H(ni,l)*drmuhl)*muhl*ww_H(nj,l))
867 !!$ !+ stab_div*mode/ray*(ww_H(ni,l)/ray+dw_H(1,ni,l,m))*ww_H(nj,l))
868 !!$ !
869 !!$
870 !!$ TH(8,ni,nj) = TH(8,ni,nj)+ rj_H(l,m) * ray* ( &
871 !!$ - mode/ray**2 * ww_H(ni,l)*(ww_H(nj,l)+ray*dw_H(1,nj,l,m))/sigma(m) &
872 !!$ !DIVERGENCE, June 8 2008
873 !!$ - c_div*mode/ray*(muhl*(ww_H(ni,l)/ray+dw_H(1,ni,l,m)) &
874 !!$ + ww_H(ni,l)*drmuhl)*muhl*ww_H(nj,l))
875 !!$ !-stab_div*mode/ray*(ww_H(ni,l)/ray+dw_H(1,ni,l,m))*ww_H(nj,l))
876 !!$ !
877 !!$
878 !!$ TH(3,ni,nj) = TH(3,ni,nj)+ rj_H(l,m) * ray* ( &
879 !!$ - dw_H(2,ni,l,m)*dw_H(1,nj,l,m)/sigma(m) &
880 !!$ !DIVERGENCE, June 8 2008
881 !!$ + c_div*(muhl*(ww_H(ni,l)/ray+dw_H(1,ni,l,m)) + ww_H(ni,l)*drmuhl)*&
882 !!$ (muhl*dw_H(2,nj,l,m) + ww_H(nj,l)*dzmuhl))
883 !!$ !+ stab_div*(ww_H(ni,l)/ray+dw_H(1,ni,l,m))*dw_H(2,nj,l,m))
884 !!$ !
885 !!$
886 !!$ TH(4,ni,nj) = TH(4,ni,nj) + rj_H(l,m) * ray* ( &
887 !!$ c_mu_H*ww_H(ni,l)*ww_H(nj,l) &
888 !!$ + (dw_H(2,ni,l,m)*dw_H(2,nj,l,m) &
889 !!$ + 1/ray**2 *(ww_H(ni,l)+ray*dw_H(1,ni,l,m))*(ww_H(nj,l)&
890 !!$ +ray*dw_H(1,nj,l,m)))/sigma(m) &
891 !!$ !DIVERGENCE, June 8 2008
892 !!$ +c_div*muhl**2*mode**2/ray**2*ww_H(ni,l)*ww_H(nj,l))
893 !!$ !+stab_div*mode**2/ray**2*ww_H(ni,l)*ww_H(nj,l))
894 !!$ !
895 !!$
896 !!$ TH(5,ni,nj) = TH(5,ni,nj) + rj_H(l,m) * ray* (&
897 !!$ + mode/ray*dw_H(2,ni,l,m)*ww_H(nj,l)/sigma(m) &
898 !!$ !DIVERGENCE, June 8 2008
899 !!$ +c_div*mode/ray*muhl*ww_H(ni,l)*(muhl*dw_H(2,nj,l,m) + ww_H(nj,l)*dzmuhl))
900 !!$ !+stab_div*mode/ray*ww_H(ni,l)*dw_H(2,nj,l,m))
901 !!$ !
902 !!$
903 !!$ TH(9,ni,nj) = TH(9,ni,nj) + rj_H(l,m) * ray* (&
904 !!$ - mode/ray*dw_H(2,ni,l,m)*ww_H(nj,l)/sigma(m) &
905 !!$ !DIVERGENCE, June 8 2008
906 !!$ - c_div*mode/ray*muhl*ww_H(ni,l)*(muhl*dw_H(2,nj,l,m) + ww_H(nj,l)*dzmuhl))
907 !!$ !- stab_div*mode/ray*ww_H(ni,l)*dw_H(2,nj,l,m))
908 !!$ !
909 !!$
910 !!$ TH(6,ni,nj) = TH(6,ni,nj) + rj_H(l,m) * ray* ( &
911 !!$ c_mu_H*ww_H(ni,l)*ww_H(nj,l) &
912 !!$ + (mode**2/ray**2*ww_H(ni,l)*ww_H(nj,l) + dw_H(1,ni,l,m)*dw_H(1,nj,l,m))/sigma(m) &
913 !!$ !DIVERGENCE, June 8 2008
914 !!$ + c_div*(muhl*dw_H(2,ni,l,m) + ww_H(ni,l)*dzmuhl) &
915 !!$ *(muhl*dw_H(2,nj,l,m) + ww_H(nj,l)*dzmuhl))
916 !!$ !+ stab_div*dw_H(2,ni,l,m)*dw_H(2,nj,l,m))
917 !!$ !
918  ENDDO
919  END DO
920 
921  END DO
922 
923  mat_loc1 = 0.d0
924  mat_loc2 = 0.d0
925  DO ki= 1, 3
926  DO ni = 1, n_wh
927  i = h_mesh%jj(ni, m)
928  ib = la_h%loc_to_glob(ki,i)
929  ix = (ki-1)*n_wh+ni
930  idxn(ix) = ib - 1
931  DO kj = 1, 3
932  DO nj = 1, n_wh
933  j = h_mesh%jj(nj, m)
934  jb = la_h%loc_to_glob(kj,j)
935  jx = (kj-1)*n_wh+nj
936  jdxn(jx) = jb - 1
937 
938  IF ((ki == 1) .AND. (kj == 1)) THEN
939  mat_loc1(ix,jx) = th(1,ni,nj)
940  mat_loc2(ix,jx) = th(1,ni,nj)
941  ELSEIF ((ki == 1) .AND. (kj == 2)) THEN
942  mat_loc1(ix,jx) = th(2,ni,nj)
943  mat_loc2(ix,jx) = th(8,ni,nj)
944  ELSEIF ((ki == 2) .AND. (kj == 1)) THEN
945  mat_loc1(ix,jx) = th(2,nj,ni)
946  mat_loc2(ix,jx) = th(8,nj,ni)
947  ELSEIF ((ki == 1) .AND. (kj == 3)) THEN
948  mat_loc1(ix,jx) = th(3,ni,nj)
949  mat_loc2(ix,jx) = th(3,ni,nj)
950  ELSEIF ((ki == 3) .AND. (kj == 1)) THEN
951  mat_loc1(ix,jx) = th(3,nj,ni)
952  mat_loc2(ix,jx) = th(3,nj,ni)
953  ELSEIF ((ki == 2) .AND. (kj == 2)) THEN
954  mat_loc1(ix,jx) = th(4,ni,nj)
955  mat_loc2(ix,jx) = th(4,ni,nj)
956  ELSEIF ((ki == 2) .AND. (kj == 3)) THEN
957  mat_loc1(ix,jx) = th(5,ni,nj)
958  mat_loc2(ix,jx) = th(9,ni,nj)
959  ELSEIF ((ki == 3) .AND. (kj == 2)) THEN
960  mat_loc1(ix,jx) = th(5,nj,ni)
961  mat_loc2(ix,jx) = th(9,nj,ni)
962  ELSEIF ((ki == 3) .AND. (kj == 3)) THEN
963  mat_loc1(ix,jx) = th(6,ni,nj)
964  mat_loc2(ix,jx) = th(6,ni,nj)
965  ENDIF
966 
967  END DO
968  END DO
969  END DO
970  END DO
971  CALL matsetvalues(h_p_phi_mat1, 3*n_wh, idxn(1:3*n_wh), 3*n_wh, jdxn(1:3*n_wh), &
972  mat_loc1(1:3*n_wh,1:3*n_wh), add_values, ierr)
973  CALL matsetvalues(h_p_phi_mat2, 3*n_wh, idxn(1:3*n_wh), 3*n_wh, jdxn(1:3*n_wh), &
974  mat_loc2(1:3*n_wh,1:3*n_wh), add_values, ierr)
975  END DO
976 
977  ! Block on Pmag
978  DO m = 1, pmag_mesh%me
979  hloc = stab_div*sqrt(sum(pmag_mesh%gauss%rj(:,m)))**(2*(1-alpha))
980  tpmag = 0.d0
981  DO l = 1, pmag_mesh%gauss%l_G
982  !Normalization
983  muhl = 1 ! SUM(mu_H_field(H_mesh%jj(1:3,m))*pmag_mesh%gauss%ww(1:3,l))
984  !ATTENTION ATTENTION: above line should be replaced by the next one
985  !JLG DCQ, July 17, 2013 (this should be the proper normalization)
986  !muhl = SUM(mu_H_field(H_mesh%jj(1:3,m))*pmag_mesh%gauss%ww(1:3,l))
987  !Normalization
988  !===Compute radius of Gauss point
989  ray = 0
990  DO ni = 1, pmag_mesh%gauss%n_w
991  i = pmag_mesh%jj(ni,m)
992  ray = ray + pmag_mesh%rr(1,i)*pmag_mesh%gauss%ww(ni,l)
993  END DO
994  viscolm = hloc*muhl*pmag_mesh%gauss%rj(l,m)
995  DO nj = 1, pmag_mesh%gauss%n_w
996  j = pmag_mesh%jj(nj, m)
997  DO ni = 1, pmag_mesh%gauss%n_w
998  i = pmag_mesh%jj(ni, m)
999  !grad(u).grad(v) en r et z
1000  xij = 0.d0
1001  DO k = 1, 2
1002  xij = xij + pmag_mesh%gauss%dw(k,nj,l,m) * pmag_mesh%gauss%dw(k,ni,l,m)
1003  END DO
1004  !blocs diagonaux
1005  tpmag(ni,nj) = tpmag(ni,nj) + ray * viscolm* xij &
1006  + viscolm*mode**2*pmag_mesh%gauss%ww(ni,l)*pmag_mesh%gauss%ww(nj,l)/ray
1007  ENDDO
1008  ENDDO
1009  ENDDO
1010 
1011  DO ni = 1, pmag_mesh%gauss%n_w
1012  i = pmag_mesh%jj(ni, m)
1013  ib = la_pmag%loc_to_glob(1,i)
1014  idxn(ni) = ib - 1
1015  DO nj = 1, pmag_mesh%gauss%n_w
1016  j = pmag_mesh%jj(nj, m)
1017  jb = la_pmag%loc_to_glob(1,j)
1018  jdxn(nj) = jb - 1
1019  END DO
1020  END DO
1021  CALL matsetvalues(h_p_phi_mat1, n_wpmag, idxn(1:n_wpmag), n_wpmag, jdxn(1:n_wpmag), &
1022  tpmag(1:n_wpmag,1:n_wpmag), add_values, ierr)
1023  CALL matsetvalues(h_p_phi_mat2, n_wpmag, idxn(1:n_wpmag), n_wpmag, jdxn(1:n_wpmag), &
1024  tpmag(1:n_wpmag,1:n_wpmag), add_values, ierr)
1025  ENDDO
1026  ! End Block on PmagxPmag
1027 
1028  ! Block on PmagxH and HxPmag
1029  DO m = 1, pmag_mesh%me
1030  IF (h_mesh%gauss%n_w==3) THEN
1031  dwp=h_mesh%gauss%dw(:,:,:,m)
1032  wwp=h_mesh%gauss%ww
1033  ELSE
1034  dwp(:,1,:) = h_mesh%gauss%dw(:,1,:,m) + 0.5d0*(h_mesh%gauss%dw(:,5,:,m)+h_mesh%gauss%dw(:,6,:,m))
1035  dwp(:,2,:) = h_mesh%gauss%dw(:,2,:,m) + 0.5d0*(h_mesh%gauss%dw(:,6,:,m)+h_mesh%gauss%dw(:,4,:,m))
1036  dwp(:,3,:) = h_mesh%gauss%dw(:,3,:,m) + 0.5d0*(h_mesh%gauss%dw(:,4,:,m)+h_mesh%gauss%dw(:,5,:,m))
1037  wwp(1,:) = h_mesh%gauss%ww(1,:) + 0.5d0*(h_mesh%gauss%ww(5,:)+h_mesh%gauss%ww(6,:))
1038  wwp(2,:) = h_mesh%gauss%ww(2,:) + 0.5d0*(h_mesh%gauss%ww(6,:)+h_mesh%gauss%ww(4,:))
1039  wwp(3,:) = h_mesh%gauss%ww(3,:) + 0.5d0*(h_mesh%gauss%ww(4,:)+h_mesh%gauss%ww(5,:))
1040  END IF
1041 
1042  thpmag = 0.d0
1043  DO l = 1, h_mesh%gauss%l_G
1044  ray = 0.d0
1045  DO ni = 1, h_mesh%gauss%n_w
1046  i = h_mesh%jj(ni,m)
1047  ray = ray + h_mesh%rr(1,i)*h_mesh%gauss%ww(ni,l)
1048  END DO
1049  muhl = stab_div*ray*h_mesh%gauss%rj(l,m)*sum(mu_h_field(h_mesh%jj(:,m))*h_mesh%gauss%ww(:,l))
1050  DO nj = 1, pmag_mesh%gauss%n_w
1051  j = pmag_mesh%jj(nj, m)
1052  DO ni = 1, h_mesh%gauss%n_w
1053  i = h_mesh%jj(ni, m)
1054  thpmag(1,ni,nj) = thpmag(1,ni,nj) + muhl*dwp(1,nj,l)*h_mesh%gauss%ww(ni,l)
1055  thpmag(2,ni,nj) = thpmag(2,ni,nj) - muhl*mode*wwp(nj,l)*h_mesh%gauss%ww(ni,l)/ray
1056  thpmag(3,ni,nj) = thpmag(3,ni,nj) + muhl*dwp(2,nj,l)*h_mesh%gauss%ww(ni,l)
1057  END DO
1058  END DO
1059  END DO
1060 
1061  mat_loc1 = 0.d0
1062  mat_loc2 = 0.d0
1063  idxn = 0
1064  jdxn = 0
1065  DO ni = 1, n_wh
1066  i = h_mesh%jj(ni, m)
1067  DO k = 1, 3
1068  IF (k==2) THEN
1069  eps=-1
1070  ELSE
1071  eps=1
1072  END IF
1073  ib = la_h%loc_to_glob(k,i)
1074  ix = (k-1)*n_wh + ni
1075  idxn(ix) = ib - 1
1076  DO nj = 1, n_wpmag
1077  j = pmag_mesh%jj(nj, m)
1078  jb = la_pmag%loc_to_glob(1,j)
1079  jx = nj
1080  jdxn(jx) = jb - 1
1081  mat_loc1(ix,jx) = thpmag(k,ni,nj)
1082  mat_loc2(ix,jx) = eps*thpmag(k,ni,nj)
1083  END DO
1084  END DO
1085  END DO
1086 
1087  CALL matsetvalues(h_p_phi_mat1, 3*n_wh, idxn(1:3*n_wh), n_wpmag, jdxn(1:n_wpmag), &
1088  mat_loc1(1:3*n_wh,1:n_wpmag), add_values, ierr)
1089  CALL matsetvalues(h_p_phi_mat2, 3*n_wh, idxn(1:3*n_wh), n_wpmag, jdxn(1:n_wpmag), &
1090  mat_loc2(1:3*n_wh,1:n_wpmag), add_values, ierr)
1091 
1092  mat_loc1 = 0.d0
1093  mat_loc2 = 0.d0
1094  DO ni = 1, n_wpmag
1095  i = pmag_mesh%jj(ni, m)
1096  ib = la_pmag%loc_to_glob(1,i)
1097  ix = ni!+3*n_wH
1098  idxn(ix) = ib - 1
1099  DO k = 1, 3
1100  IF (k==2) THEN
1101  eps=-1
1102  ELSE
1103  eps=1
1104  END IF
1105  DO nj = 1, n_wh
1106  j = h_mesh%jj(nj, m)
1107  jb = la_h%loc_to_glob(k,j)
1108  jx = (k-1)*n_wh + nj
1109  jdxn(jx) = jb - 1
1110  mat_loc1(ix,jx) = - thpmag(k,nj,ni)
1111  mat_loc2(ix,jx) = - eps*thpmag(k,nj,ni)
1112  END DO
1113  END DO
1114  END DO
1115  CALL matsetvalues(h_p_phi_mat1, n_wpmag, idxn(1:n_wpmag), 3*n_wh, jdxn(1:3*n_wh), &
1116  mat_loc1(1:n_wpmag,1:3*n_wh), add_values, ierr)
1117  CALL matsetvalues(h_p_phi_mat2, n_wpmag, idxn(1:n_wpmag), 3*n_wh, jdxn(1:3*n_wh), &
1118  mat_loc2(1:n_wpmag,1:3*n_wh), add_values, ierr)
1119  END DO
1120  ! End Block on PmagxH and HxPmag
1121 
1122  !==Block on phi
1123  DO m = 1,phi_mesh%me
1124 
1125  tphi = 0.d0
1126 
1127  DO l = 1, phi_mesh%gauss%l_G
1128 
1129  !===Compute radius of Gauss point
1130  ray = 0
1131  DO ni = 1, phi_mesh%gauss%n_w; i = phi_mesh%jj(ni,m)
1132  ray = ray + phi_mesh%rr(1,i)*ww_phi(ni,l)
1133  END DO
1134 
1135  DO ni = 1, phi_mesh%gauss%n_w
1136  DO nj = 1, phi_mesh%gauss%n_w
1137 
1138  !mu_phi * <Grad bi, Grad bj>
1139  !JLG, FL May 28, 2009
1140  !On ajoute le laplacien de phi.
1141  !TPhi(ni,nj) = TPhi(ni,nj) + rj_phi(l,m) * ray* (c_mu_phi) &
1142  ! *(dw_phi(1,ni,l,m)*dw_phi(1,nj,l,m)+dw_phi(2,ni,l,m)*dw_phi(2,nj,l,m) &
1143  ! +mode**2/ray**2*ww_phi(ni,l)*ww_phi(nj,l))
1144  tphi(ni,nj) = tphi(ni,nj) + rj_phi(l,m) * ray* (c_mass+c_lap)*mu_phi &
1145  *(dw_phi(1,ni,l,m)*dw_phi(1,nj,l,m)+dw_phi(2,ni,l,m)*dw_phi(2,nj,l,m) &
1146  +mode**2/ray**2*ww_phi(ni,l)*ww_phi(nj,l))
1147  !JLG, FL May 28, 2009
1148  ENDDO
1149  END DO
1150  END DO
1151 
1152  !TEST
1153  !TPhi = 0.d0
1154  !TEST
1155 
1156  DO ni = 1, phi_mesh%gauss%n_w
1157  i = phi_mesh%jj(ni, m)
1158  ib = la_phi%loc_to_glob(1,i)
1159  idxn(ni) = ib - 1
1160  DO nj = 1, phi_mesh%gauss%n_w
1161  j = phi_mesh%jj(nj, m)
1162  jb = la_phi%loc_to_glob(1,j)
1163  jdxn(nj) = jb - 1
1164  END DO
1165  END DO
1166  CALL matsetvalues(h_p_phi_mat1, n_wphi, idxn(1:n_wphi), n_wphi, jdxn(1:n_wphi), &
1167  tphi(1:n_wphi,1:n_wphi), add_values, ierr)
1168  CALL matsetvalues(h_p_phi_mat2, n_wphi, idxn(1:n_wphi), n_wphi, jdxn(1:n_wphi), &
1169  tphi(1:n_wphi,1:n_wphi), add_values, ierr)
1170  END DO
1171 
1172  !*********************************************************************************
1173  !--------------------TERMS on interface_H_phi SIGMA-------------------------------
1174  !**********************************************************************************
1175 
1176  !WRITE(*,*) 'Assembling interface_H_phi '
1177  CALL gauss(phi_mesh)
1178  n_ws1 = h_mesh%gauss%n_ws
1179  n_ws2 = phi_mesh%gauss%n_ws
1180  n_w1 = h_mesh%gauss%n_w
1181  n_w2 = phi_mesh%gauss%n_w
1182 
1183  IF (h_mesh%gauss%n_ws == n_ws) THEN
1184 
1185  DO ms = 1, interface_h_phi%mes
1186 
1187  ms2 = interface_h_phi%mesh2(ms)
1188  m2 = phi_mesh%neighs(ms2)
1189  ms1 = interface_h_phi%mesh1(ms)
1190  m1 = h_mesh%neighs(ms1)
1191 
1192  ref = 1.d-8+sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - h_mesh%rr(:,h_mesh%jjs(2,ms1)))**2)
1193  diff = sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - phi_mesh%rr(:,phi_mesh%jjs(1,ms2)))**2)
1194  IF (diff/ref .LT. 1.d-10) THEN ! 1 = 1
1195  w_cs = wws
1196  ELSE ! 1 = 2
1197  DO ls = 1, l_gs
1198  w_cs(1,ls)= wws(2,ls)
1199  w_cs(2,ls)= wws(1,ls)
1200  w_cs(3,ls)= wws(3,ls)
1201  WRITE(*,*) ' Ouaps! oder of shape functions changed?'
1202  END DO
1203  END IF
1204 
1205  DO ls = 1, l_gs
1206  gauss2(1,ls) = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))*phi_mesh%gauss%wws(:,ls))
1207  gauss2(2,ls) = sum(phi_mesh%rr(2,phi_mesh%jjs(:,ms2))*phi_mesh%gauss%wws(:,ls))
1208  gauss1(1,ls) = sum( h_mesh%rr(1, h_mesh%jjs(:,ms1))* h_mesh%gauss%wws(:,ls))
1209  gauss1(2,ls) = sum( h_mesh%rr(2, h_mesh%jjs(:,ms1))* h_mesh%gauss%wws(:,ls))
1210  END DO
1211 
1212  DO ls2 = 1, l_gs
1213  ref = sqrt(1.d-8+sum(gauss2(:,ls2)**2))
1214  mark = .false.
1215  DO ls1 = 1, l_gs
1216  diff = sqrt(sum((gauss2(:,ls2)-gauss1(:,ls1))**2))
1217  IF (diff .LT. 1.d-10) THEN
1218  dw_cs(:,:,ls2,ms1) = h_mesh%gauss%dw_s(:,:,ls1,ms1)
1219  mark = .true.
1220  EXIT
1221  END IF
1222  END DO
1223  IF (.NOT.mark) WRITE(*,*) ' BUG '
1224  END DO
1225 
1226  END DO
1227 
1228  ELSE
1229  DO ms = 1, interface_h_phi%mes
1230 
1231  ms2 = interface_h_phi%mesh2(ms)
1232  m2 = phi_mesh%neighs(ms2)
1233  ms1 = interface_h_phi%mesh1(ms)
1234  m1 = h_mesh%neighs(ms1)
1235 
1236  ref = 1.d-8+sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - h_mesh%rr(:,h_mesh%jjs(2,ms1)))**2)
1237  diff = sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - phi_mesh%rr(:,phi_mesh%jjs(1,ms2)))**2)
1238  IF (diff/ref .LT. 1.d-10) THEN ! 1 = 1
1239  DO ls = 1, l_gs
1240  w_cs(1,ls)= wws(1,ls)+0.5*wws(3,ls)
1241  w_cs(2,ls)= wws(2,ls)+0.5*wws(3,ls)
1242  w_cs(3,ls)= 0
1243  END DO
1244  ELSE ! 1 = 2
1245  DO ls = 1, l_gs
1246  w_cs(1,ls)= wws(2,ls)+0.5*wws(3,ls)
1247  w_cs(2,ls)= wws(1,ls)+0.5*wws(3,ls)
1248  w_cs(3,ls)= 0
1249  WRITE(*,*) ' Ouaps! oder of shape functions changed?'
1250  END DO
1251  END IF
1252 
1253  DO ls = 1, l_gs
1254  dw_cs(1,:,ls,ms1) = h_mesh%gauss%dw(1,:,1,m1)
1255  dw_cs(2,:,ls,ms1) = h_mesh%gauss%dw(2,:,1,m1)
1256  END DO
1257 
1258  END DO
1259  END IF
1260 
1261  error = 0
1262  DO ms = 1, interface_h_phi%mes
1263 
1264  ms2 = interface_h_phi%mesh2(ms)
1265  ms1 = interface_h_phi%mesh1(ms)
1266  m2 = phi_mesh%neighs(ms2)
1267  m1 = h_mesh%neighs(ms1)
1268  mu_h = sum(mu_h_field(h_mesh%jj(:,m1)))/h_mesh%gauss%n_w
1269  !JLG, FL, May, 28, 2009
1270  hm1 = stab_colle_h_phi/sum(rjs(:,ms2))
1271  !hm1 = stab_colle_H_phi*(((mu_phi+mu_H)/mu_H)/SUM(rjs(:,ms2)))
1272  !JLG, FL, May, 28, 2009
1273 
1274  !====================================================================================
1275  !------------------------------------TERMES SUR LE BLOC H----------------------------
1276  !====================================================================================
1277 
1278  !-------------------------------hm1 (bi x ni) . (bj x nj)----------------------------
1279  !====================================================================================
1280 
1281  hsij = 0.d0
1282  DO ls = 1, l_gs
1283  !===Compute radius of Gauss point
1284  ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1285  x = hm1*rjs(ls,ms2)*ray
1286 
1287  DO ni = 1, n_ws1
1288  DO nj = 1, n_ws1
1289  y = x * w_cs(ni,ls)*w_cs(nj,ls)
1290  hsij(1,ni,nj) = hsij(1,ni,nj) + y*(rnorms(2,ls,ms2)**2)
1291  hsij(4,ni,nj) = hsij(4,ni,nj) - y*rnorms(1,ls,ms2)*rnorms(2,ls,ms2)
1292  hsij(5,ni,nj) = hsij(5,ni,nj) + y
1293  hsij(6,ni,nj) = hsij(6,ni,nj) + y*(rnorms(1,ls,ms2)**2)
1294  ENDDO
1295  ENDDO
1296 
1297  ENDDO
1298 
1299 
1300  !TEST
1301  !Hsij = 0.d0
1302  !Hsij = Hsij / hm1
1303  !TEST
1304  mat_loc1 = 0.d0
1305  mat_loc2 = 0.d0
1306  DO ki= 1, 3
1307  DO ni = 1, n_ws1
1308  i = interface_h_phi%jjs1(ni,ms)
1309  ib = la_h%loc_to_glob(ki,i)
1310  ix = (ki-1)*n_ws1+ni
1311  idxn(ix) = ib - 1
1312  DO kj = 1, 3
1313  DO nj = 1, n_ws1
1314  j = interface_h_phi%jjs1(nj,ms)
1315  jb = la_h%loc_to_glob(kj,j)
1316  jx = (kj-1)*n_ws1+nj
1317  jdxn(jx) = jb - 1
1318  IF ((ki == 1) .AND. (kj == 1)) THEN
1319  mat_loc1(ix,jx) = hsij(1,ni,nj)
1320  mat_loc2(ix,jx) = hsij(1,ni,nj)
1321  ELSEIF ((ki == 1) .AND. (kj == 3)) THEN
1322  mat_loc1(ix,jx) = hsij(4,ni,nj)
1323  mat_loc2(ix,jx) = hsij(4,ni,nj)
1324  ELSEIF ((ki == 3) .AND. (kj == 1)) THEN
1325  mat_loc1(ix,jx) = hsij(4,nj,ni)
1326  mat_loc2(ix,jx) = hsij(4,nj,ni)
1327  ELSEIF ((ki == 2) .AND. (kj == 2)) THEN
1328  mat_loc1(ix,jx) = hsij(5,ni,nj)
1329  mat_loc2(ix,jx) = hsij(5,ni,nj)
1330  ELSEIF ((ki == 3) .AND. (kj == 3)) THEN
1331  mat_loc1(ix,jx) = hsij(6,ni,nj)
1332  mat_loc2(ix,jx) = hsij(6,ni,nj)
1333  ENDIF
1334  END DO
1335  END DO
1336  END DO
1337  END DO
1338  CALL matsetvalues(h_p_phi_mat1, 3*n_ws1, idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
1339  mat_loc1(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
1340  CALL matsetvalues(h_p_phi_mat2, 3*n_ws1, idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
1341  mat_loc2(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
1342 
1343  !====================================================================================
1344  !------------------------(1/sigma) (Rot bj) . (bi x ni)------------------------------
1345  !====================================================================================
1346 
1347  hsij = 0.d0
1348  DO ls = 1, phi_mesh%gauss%l_Gs
1349 
1350  !===Compute radius of Gauss point
1351  ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1352 !!$ x = rjs(ls,ms2)*ray/sigma(m1)
1353 ! TEST DEBUG
1354  x = rjs(ls,ms2)*ray/sum(sigma_np(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
1355 ! TEST DEBUG
1356  !terme sans derivees
1357  DO ni = 1,n_ws1
1358  DO nj = 1, n_ws1
1359  y = x*w_cs(ni,ls)*w_cs(nj,ls)
1360  hsij(2,ni,nj) = hsij(2,ni,nj) + y * (-mode/ray)*(-rnorms(1,ls,ms2))
1361  hsij(3,ni,nj) = hsij(3,ni,nj) + y * mode/ray *(-rnorms(1,ls,ms2))
1362  hsij(5,ni,nj) = hsij(5,ni,nj) + y * (-1/ray) *(-rnorms(1,ls,ms2))
1363  hsij(8,ni,nj) = hsij(8,ni,nj) + y * (-mode/ray)*(-rnorms(2,ls,ms2))
1364  hsij(9,ni,nj) = hsij(9,ni,nj) + y * mode/ray *(-rnorms(2,ls,ms2))
1365  ENDDO
1366  ENDDO
1367 
1368  ENDDO
1369 
1370  !TEST
1371  !Hsij = 0.d0
1372  !TEST
1373 
1374  mat_loc1 = 0.d0
1375  mat_loc2 = 0.d0
1376  DO ki= 1, 3
1377  DO ni = 1, n_ws1
1378  i = interface_h_phi%jjs1(ni,ms)
1379  ib = la_h%loc_to_glob(ki,i)
1380  ix = (ki-1)*n_ws1 + ni
1381  idxn(ix) = ib - 1
1382  DO kj = 1, 3
1383  DO nj = 1, n_ws1
1384  j = interface_h_phi%jjs1(nj,ms)
1385  jb = la_h%loc_to_glob(kj,j)
1386  jx = (kj-1)*n_ws1 + nj
1387  jdxn(jx) = jb - 1
1388  IF ( (ki == 2) .AND. (kj == 1)) THEN
1389  mat_loc1(ix,jx) = hsij(2,ni,nj)
1390  mat_loc2(ix,jx) = hsij(3,ni,nj)
1391  ELSEIF ((ki == 2) .AND. (kj == 2)) THEN
1392  mat_loc1(ix,jx) = hsij(5,ni,nj)
1393  mat_loc2(ix,jx) = hsij(5,ni,nj)
1394  ELSEIF ( (ki == 2) .AND. (kj == 3)) THEN
1395  mat_loc1(ix,jx) = hsij(8,ni,nj)
1396  mat_loc2(ix,jx) = hsij(9,ni,nj)
1397  ENDIF
1398  END DO
1399  END DO
1400  END DO
1401  END DO
1402  CALL matsetvalues(h_p_phi_mat1, 3*n_ws1, idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
1403  mat_loc1(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
1404  CALL matsetvalues(h_p_phi_mat2, 3*n_ws1, idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
1405  mat_loc2(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
1406 
1407  !Feb 2 2007
1408  mat_loc1 = 0.d0
1409  mat_loc2 = 0.d0
1410  hsij=c_sym*hsij !SYM
1411  DO ki= 1, 3
1412  DO ni = 1, n_ws1
1413  i = interface_h_phi%jjs1(ni,ms)
1414  ib = la_h%loc_to_glob(ki,i)
1415  ix = (ki-1)*n_ws1 + ni
1416  idxn(ix) = ib - 1
1417  DO kj = 1, 3
1418  DO nj = 1, n_ws1
1419  j = interface_h_phi%jjs1(nj,ms)
1420  jb = la_h%loc_to_glob(kj,j)
1421  jx = (kj-1)*n_ws1 + nj
1422  jdxn(jx) = jb - 1
1423  IF ( (kj == 2) .AND. (ki == 1)) THEN
1424  mat_loc1(ix,jx) = hsij(2,nj,ni)
1425  mat_loc2(ix,jx) = hsij(3,nj,ni)
1426  ELSEIF ((kj == 2) .AND. (ki == 2)) THEN
1427  mat_loc1(ix,jx) = hsij(5,nj,ni)
1428  mat_loc2(ix,jx) = hsij(5,nj,ni)
1429  ELSEIF ( (kj == 2) .AND. (ki == 3)) THEN
1430  mat_loc1(ix,jx) = hsij(8,nj,ni)
1431  mat_loc2(ix,jx) = hsij(9,nj,ni)
1432  ENDIF
1433  END DO
1434  END DO
1435  END DO
1436  END DO
1437  !feb 2 2007
1438  CALL matsetvalues(h_p_phi_mat1, 3*n_ws1, idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
1439  mat_loc1(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
1440  CALL matsetvalues(h_p_phi_mat2, 3*n_ws1, idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
1441  mat_loc2(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
1442 
1443  hsij = 0.d0
1444  DO ls = 1, phi_mesh%gauss%l_Gs
1445 
1446  !===Compute radius of Gauss point
1447  ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1448 !!$ x = rjs(ls,ms2)*ray /sigma(m1)
1449 ! TEST DEBUG
1450  x = rjs(ls,ms2)*ray/sum(sigma_np(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
1451 ! TEST DEBUG
1452  !termes avec derivees
1453  DO ni = 1,n_ws1
1454  y = x*w_cs(ni,ls)
1455  DO nj = 1, n_w1
1456  hsij(1,ni,nj) = hsij(1,ni,nj) + y*(-dw_cs(2,nj,ls,ms1))*(-rnorms(2,ls,ms2))
1457  hsij(4,ni,nj) = hsij(4,ni,nj) + y* dw_cs(1,nj,ls,ms1) *(-rnorms(2,ls,ms2))
1458  hsij(5,ni,nj) = hsij(5,ni,nj) + &
1459  y*(-dw_cs(2,nj,ls,ms1)*(-rnorms(2,ls,ms2))-dw_cs(1,nj,ls,ms1)*(-rnorms(1,ls,ms2)))
1460  hsij(6,ni,nj) = hsij(6,ni,nj) + y*(-dw_cs(1,nj,ls,ms1))*(-rnorms(1,ls,ms2))
1461  hsij(7,ni,nj) = hsij(7,ni,nj) + y* dw_cs(2,nj,ls,ms1) *(-rnorms(1,ls,ms2))
1462  ENDDO
1463  ENDDO
1464 
1465  ENDDO
1466 
1467  !TEST
1468  !Hsij = 0.d0
1469  !TEST
1470  mat_loc1 = 0.d0
1471  mat_loc2 = 0.d0
1472  DO ki= 1, 3
1473  DO ni = 1, n_ws1
1474  i = interface_h_phi%jjs1(ni,ms)
1475  ib = la_h%loc_to_glob(ki,i)
1476  ix = (ki-1)*n_ws1 + ni
1477  idxn(ix) = ib - 1
1478  DO kj = 1, 3
1479  DO nj = 1, n_w1
1480  j = h_mesh%jj(nj,m1)
1481  jb = la_h%loc_to_glob(kj,j)
1482  jx = (kj-1)*n_w1 + nj
1483  jdxn(jx) = jb - 1
1484  IF ((ki == 1) .AND. (kj == 1)) THEN
1485  mat_loc1(ix,jx) = hsij(1,ni,nj)
1486  mat_loc2(ix,jx) = hsij(1,ni,nj)
1487  ELSEIF ((ki == 1) .AND. (kj == 3)) THEN
1488  mat_loc1(ix,jx) = hsij(4,ni,nj)
1489  mat_loc2(ix,jx) = hsij(4,ni,nj)
1490  ELSEIF ((ki == 2) .AND. (kj == 2)) THEN
1491  mat_loc1(ix,jx) = hsij(5,ni,nj)
1492  mat_loc2(ix,jx) = hsij(5,ni,nj)
1493  ELSEIF ((ki == 3) .AND. (kj == 3)) THEN
1494  mat_loc1(ix,jx) = hsij(6,ni,nj)
1495  mat_loc2(ix,jx) = hsij(6,ni,nj)
1496  ELSEIF ((ki == 3) .AND. (kj == 1)) THEN
1497  mat_loc1(ix,jx) = hsij(7,ni,nj)
1498  mat_loc2(ix,jx) = hsij(7,ni,nj)
1499  ENDIF
1500  END DO
1501  END DO
1502  END DO
1503  END DO
1504 
1505  CALL matsetvalues(h_p_phi_mat1, 3*n_ws1, idxn(1:3*n_ws1), 3*n_w1, jdxn(1:3*n_w1), &
1506  mat_loc1(1:3*n_ws1,1:3*n_w1), add_values, ierr)
1507  CALL matsetvalues(h_p_phi_mat2, 3*n_ws1, idxn(1:3*n_ws1), 3*n_w1, jdxn(1:3*n_w1), &
1508  mat_loc2(1:3*n_ws1,1:3*n_w1), add_values, ierr)
1509 
1510  !Feb 2 2007
1511  mat_loc1 = 0.d0
1512  mat_loc2 = 0.d0
1513  hsij=c_sym*hsij !SYM
1514  DO ki = 1, 3
1515  DO ni = 1, n_w1
1516  i = h_mesh%jj(ni,m1)
1517  ib = la_h%loc_to_glob(ki,i)
1518  ix = (ki-1)*n_w1 + ni
1519  idxn(ix) = ib - 1
1520  DO kj= 1, 3
1521  DO nj = 1, n_ws1
1522  j = interface_h_phi%jjs1(nj,ms)
1523  jb = la_h%loc_to_glob(kj,j)
1524  jx = (kj-1)*n_ws1 + nj
1525  jdxn(jx) = jb - 1
1526  IF ((kj == 1) .AND. (ki == 1)) THEN
1527  mat_loc1(ix,jx) = hsij(1,nj,ni)
1528  mat_loc2(ix,jx) = hsij(1,nj,ni)
1529  ELSEIF ((kj == 1) .AND. (ki == 3)) THEN
1530  mat_loc1(ix,jx) = hsij(4,nj,ni)
1531  mat_loc2(ix,jx) = hsij(4,nj,ni)
1532  ELSEIF ((kj == 2) .AND. (ki == 2)) THEN
1533  mat_loc1(ix,jx) = hsij(5,nj,ni)
1534  mat_loc2(ix,jx) = hsij(5,nj,ni)
1535  ELSEIF ((kj == 3) .AND. (ki == 3)) THEN
1536  mat_loc1(ix,jx) = hsij(6,nj,ni)
1537  mat_loc2(ix,jx) = hsij(6,nj,ni)
1538  ELSEIF ((kj == 3) .AND. (ki == 1)) THEN
1539  mat_loc1(ix,jx) = hsij(7,nj,ni)
1540  mat_loc2(ix,jx) = hsij(7,nj,ni)
1541  ENDIF
1542  END DO
1543  END DO
1544  END DO
1545  END DO
1546  CALL matsetvalues(h_p_phi_mat1, 3*n_w1, idxn(1:3*n_w1), 3*n_ws1, jdxn(1:3*n_ws1), &
1547  mat_loc1(1:3*n_w1,1:3*n_ws1), add_values, ierr)
1548  CALL matsetvalues(h_p_phi_mat2, 3*n_w1, idxn(1:3*n_w1), 3*n_ws1, jdxn(1:3*n_ws1), &
1549  mat_loc2(1:3*n_w1,1:3*n_ws1), add_values, ierr)
1550  !Feb 2 2007
1551 
1552 
1553  !====================================================================================
1554  !------------------------------------TERMES SUR LE BLOC PHI--------------------------
1555  !====================================================================================
1556 
1557  !------------------------hm1 (Grad(phi_i) x ni).(Grad(phi_j) x nj)-------------------
1558  !====================================================================================
1559 
1560  phisij = 0.d0
1561 
1562  DO ls = 1, phi_mesh%gauss%l_Gs
1563 
1564  !===Compute radius of Gauss point
1565  ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1566  x = hm1*rjs(ls,ms2)*ray
1567 
1568  !terme sans derivee
1569  DO ni=1, n_ws2
1570  DO nj=1, n_ws2
1571  phisij(ni,nj) = phisij(ni,nj) + x*mode**2/ray**2*wws(ni,ls)*wws(nj,ls)
1572  ENDDO
1573  ENDDO
1574 
1575  ENDDO
1576 
1577  !TEST
1578  !Phisij = 0.d0
1579  !Phisij = Phisij/hm1
1580  !TEST
1581  DO ni = 1, n_ws2
1582  i = interface_h_phi%jjs2(ni,ms)
1583  ib = la_phi%loc_to_glob(1,i)
1584  idxn(ni) = ib - 1
1585  DO nj = 1, n_ws2
1586  j = interface_h_phi%jjs2(nj,ms)
1587  jb = la_phi%loc_to_glob(1,j)
1588  jdxn(nj) = jb - 1
1589  END DO
1590  END DO
1591  CALL matsetvalues(h_p_phi_mat1, n_ws2, idxn(1:n_ws2), n_ws2, jdxn(1:n_ws2), &
1592  phisij(1:n_ws2,1:n_ws2), add_values, ierr)
1593  CALL matsetvalues(h_p_phi_mat2, n_ws2, idxn(1:n_ws2), n_ws2, jdxn(1:n_ws2), &
1594  phisij(1:n_ws2,1:n_ws2), add_values, ierr)
1595 
1596  phisij = 0.d0
1597  DO ls = 1, l_gs
1598 
1599  !===Compute radius of Gauss point
1600  ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1601  x = hm1*rjs(ls,ms2)*ray
1602 
1603  !terme avec derivee
1604  DO ni = 1, n_w2
1605  DO nj = 1, n_w2
1606  phisij(ni,nj) = phisij(ni,nj) + x*( &
1607  (dw_s(2,ni,ls,ms2)*rnorms(1,ls,ms2) - dw_s(1,ni,ls,ms2)*rnorms(2,ls,ms2))* &
1608  (dw_s(2,nj,ls,ms2)*rnorms(1,ls,ms2) - dw_s(1,nj,ls,ms2)*rnorms(2,ls,ms2)))
1609  ENDDO
1610  ENDDO
1611  ENDDO
1612 
1613  !Phisij = 0.d0
1614  !Phisij = Phisij/hm1
1615  !TEST
1616 
1617  DO ni = 1, n_w2
1618  i = phi_mesh%jj(ni, m2)
1619  ib = la_phi%loc_to_glob(1,i)
1620  idxn(ni) = ib - 1
1621  DO nj = 1, n_w2
1622  j = phi_mesh%jj(nj, m2)
1623  jb = la_phi%loc_to_glob(1,j)
1624  jdxn(nj) = jb - 1
1625  END DO
1626  END DO
1627  CALL matsetvalues(h_p_phi_mat1, n_w2, idxn(1:n_w2), n_w2, jdxn(1:n_w2), &
1628  phisij(1:n_w2,1:n_w2), add_values, ierr)
1629  CALL matsetvalues(h_p_phi_mat2, n_w2, idxn(1:n_w2), n_w2, jdxn(1:n_w2), &
1630  phisij(1:n_w2,1:n_w2), add_values, ierr)
1631  !====================================================================================
1632  !------------------------------------TERMES CROISES----------------------------------
1633  !====================================================================================
1634 
1635  !====================================================================================
1636  !------------------------hm1 (bi x ni) . (Grad(phi_j) x nj)--------------------------
1637  !------------------ + hm1(Grad(phi_i) x ni).(bj x nj)---------------------------
1638  !====================================================================================
1639 
1640  sij = 0.d0
1641  DO ls = 1, l_gs
1642 
1643  !===Compute radius of Gauss point
1644  ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1645  x = hm1*rjs(ls,ms2)*ray
1646 
1647  !terme sans derivee
1648  DO ni = 1, n_ws1
1649  DO nj = 1, n_ws2
1650  sij(3,ni,nj) = sij(3,ni,nj) + x*(mode/ray)*w_cs(ni,ls)*wws(nj,ls)
1651  ENDDO
1652  ENDDO
1653  ENDDO
1654  sij(4,:,:) = -sij(3,:,:)
1655 
1656  !TEST
1657  !Sij = 0.d0
1658  !Sij = Sij /hm1
1659  !TEST
1660 
1661  ki = 2
1662  DO ni = 1, n_ws1
1663  i = interface_h_phi%jjs1(ni,ms)
1664  ib = la_h%loc_to_glob(ki,i)
1665  idxn(ni) = ib - 1
1666  DO nj = 1, n_ws2
1667  j = interface_h_phi%jjs2(nj,ms)
1668  jb = la_phi%loc_to_glob(1,j)
1669  jdxn(nj) = jb - 1
1670  END DO
1671  ENDDO
1672  CALL matsetvalues(h_p_phi_mat1, n_ws1, idxn(1:n_ws1), n_ws2, jdxn(1:n_ws2), &
1673  sij(3,1:n_ws1,1:n_ws2), add_values, ierr)
1674  CALL matsetvalues(h_p_phi_mat2, n_ws1, idxn(1:n_ws1), n_ws2, jdxn(1:n_ws2), &
1675  sij(4,1:n_ws1,1:n_ws2), add_values, ierr)
1676 
1677  !TEST SYM
1678  !Feb 2 2003
1679  !Sij = 0.d0
1680  mat_loc1 = 0.d0
1681  mat_loc2 = 0.d0
1682  kj = 2
1683  DO ni = 1, n_ws2
1684  i = interface_h_phi%jjs2(ni,ms)
1685  ib = la_phi%loc_to_glob(1,i)
1686  idxn(ni) = ib - 1
1687  DO nj = 1, n_ws1
1688  j = interface_h_phi%jjs1(nj,ms)
1689  jb = la_h%loc_to_glob(kj,j)
1690  jdxn(nj) = jb - 1
1691  mat_loc1(ni,nj) = sij(3,nj,ni)
1692  mat_loc2(ni,nj) = sij(4,nj,ni)
1693  END DO
1694  ENDDO
1695  CALL matsetvalues(h_p_phi_mat1, n_ws2, idxn(1:n_ws2), n_ws1, jdxn(1:n_ws1), &
1696  mat_loc1(1:n_ws2,1:n_ws1), add_values, ierr)
1697  CALL matsetvalues(h_p_phi_mat2, n_ws2, idxn(1:n_ws2), n_ws1, jdxn(1:n_ws1), &
1698  mat_loc2(1:n_ws2,1:n_ws1), add_values, ierr)
1699 
1700  !Feb 2 2003
1701  !TEST SYM
1702  sij = 0.d0
1703  DO ls = 1, l_gs
1704 
1705  !===Compute radius of Gauss point
1706  ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1707  x = hm1*rjs(ls,ms2)*ray
1708 
1709  !terme avec derivee
1710  DO ni = 1, n_ws1
1711  y = x * w_cs(ni,ls)
1712  DO nj = 1, n_w2
1713  sij(1,ni,nj) = sij(1,ni,nj) + &
1714  y*(-dw_s(1,nj,ls,ms2)*rnorms(2,ls,ms2)**2 + dw_s(2,nj,ls,ms2)*rnorms(1,ls,ms2)*rnorms(2,ls,ms2))
1715  sij(5,ni,nj) = sij(5,ni,nj) + &
1716  y*(-dw_s(2,nj,ls,ms2)*rnorms(1,ls,ms2)**2 + dw_s(1,nj,ls,ms2)*rnorms(1,ls,ms2)*rnorms(2,ls,ms2))
1717  ENDDO
1718  ENDDO
1719  ENDDO
1720 
1721  !TEST
1722  !Sij = 0.d0
1723  !Sij = Sij /hm1
1724  !TEST
1725  mat_loc1 = 0.d0
1726  mat_loc2 = 0.d0
1727  DO ki= 1, 3
1728  DO ni = 1, n_ws1
1729  i = interface_h_phi%jjs1(ni,ms)
1730  ib = la_h%loc_to_glob(ki,i)
1731  ix = (ki-1)*n_ws1 + ni
1732  idxn(ix) = ib - 1
1733  DO nj = 1, n_w2
1734  j = phi_mesh%jj(nj,m2)
1735  jb = la_phi%loc_to_glob(1,j)
1736  jx = nj
1737  jdxn(jx) = jb - 1
1738  IF (ki == 1) THEN
1739  mat_loc1(ix,jx) = sij(1,ni,nj)
1740  mat_loc2(ix,jx) = sij(1,ni,nj)
1741  ELSEIF (ki == 3) THEN
1742  mat_loc1(ix,jx) = sij(5,ni,nj)
1743  mat_loc2(ix,jx) = sij(5,ni,nj)
1744  END IF
1745  END DO
1746  END DO
1747  END DO
1748  CALL matsetvalues(h_p_phi_mat1, 3*n_ws1, idxn(1:3*n_ws1), n_w2, jdxn(1:n_w2), &
1749  mat_loc1(1:3*n_ws1,1:n_w2), add_values, ierr)
1750  CALL matsetvalues(h_p_phi_mat2, 3*n_ws1, idxn(1:3*n_ws1), n_w2, jdxn(1:n_w2), &
1751  mat_loc2(1:3*n_ws1,1:n_w2), add_values, ierr)
1752 
1753  !TEST SYM
1754  !Feb 2 2003
1755  !Sij = 0.d0
1756  mat_loc1 = 0.d0
1757  mat_loc2 = 0.d0
1758  DO ni = 1, n_w2
1759  i = phi_mesh%jj(ni,m2)
1760  ib = la_phi%loc_to_glob(1,i)
1761  ix = ni
1762  idxn(ix) = ib - 1
1763  DO kj=1,3
1764  DO nj = 1, n_ws1
1765  j = interface_h_phi%jjs1(nj,ms)
1766  jb = la_h%loc_to_glob(kj,j)
1767  jx = (kj-1)*n_ws1 + nj
1768  jdxn(jx) = jb - 1
1769  IF (kj == 1) THEN
1770  mat_loc1(ix,jx) = sij(1,nj,ni)
1771  mat_loc2(ix,jx) = sij(1,nj,ni)
1772  ELSEIF (kj == 3) THEN
1773  mat_loc1(ix,jx) = sij(5,nj,ni)
1774  mat_loc2(ix,jx) = sij(5,nj,ni)
1775  ENDIF
1776  END DO
1777  END DO
1778  ENDDO
1779  CALL matsetvalues(h_p_phi_mat1, n_w2, idxn(1:n_w2), 3*n_ws1, jdxn(1:3*n_ws1), &
1780  mat_loc1(1:n_w2,1:3*n_ws1), add_values, ierr)
1781  CALL matsetvalues(h_p_phi_mat2, n_w2, idxn(1:n_w2), 3*n_ws1, jdxn(1:3*n_ws1), &
1782  mat_loc2(1:n_w2,1:3*n_ws1), add_values, ierr)
1783 
1784  !TEST SYM
1785  !Feb 2 2003
1786 
1787  !====================================================================================
1788  !----------------------(1/sigma) (Rot bj).(Grad(phi_i) x ni)-------------------------
1789  !====================================================================================
1790  ! GOTO 200
1791 
1792 
1793  sij = 0.d0
1794  DO ls = 1, l_gs
1795 
1796  !===Compute radius of Gauss point
1797  ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1798 !!$ x = rjs(ls,ms2)*ray/sigma(m1)
1799 ! TEST DEBUG
1800  x = rjs(ls,ms2)*ray/sum(sigma_np(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
1801 ! TEST DEBUG
1802  !terme sans derivee
1803  DO ni = 1, n_ws2
1804  DO nj = 1, n_ws1
1805  y = x * wws(ni,ls)*w_cs(nj,ls)
1806  sij(1,ni,nj) = sij(1,ni,nj) + y*( mode/ray)**2*rnorms(1,ls,ms2)
1807  sij(3,ni,nj) = sij(3,ni,nj) + y*( mode/ray**2)*rnorms(1,ls,ms2)
1808  sij(4,ni,nj) = sij(4,ni,nj) + y*(-mode/ray**2)*rnorms(1,ls,ms2)
1809  sij(5,ni,nj) = sij(5,ni,nj) + y*( mode/ray)**2*rnorms(2,ls,ms2)
1810  ENDDO
1811  ENDDO
1812 
1813  ENDDO
1814 
1815  !TEST
1816  !Sij = 0.d0
1817  !TEST
1818  mat_loc1 = 0.d0
1819  mat_loc2 = 0.d0
1820  DO ni = 1, n_ws2
1821  i = interface_h_phi%jjs2(ni,ms)
1822  ib = la_phi%loc_to_glob(1,i)
1823  ix = ni
1824  idxn(ix) = ib - 1
1825  DO kj =1,3
1826  DO nj = 1, n_ws1
1827  j = interface_h_phi%jjs1(nj,ms)
1828  jb = la_h%loc_to_glob(kj,j)
1829  jx = (kj-1)*n_ws1 + nj
1830  jdxn(jx) = jb - 1
1831  IF (kj == 1) THEN
1832  mat_loc1(ix,jx) = sij(1,ni,nj)
1833  mat_loc2(ix,jx) = sij(1,ni,nj)
1834  ELSEIF (kj == 2) THEN
1835  mat_loc1(ix,jx) = sij(3,ni,nj)
1836  mat_loc2(ix,jx) = sij(4,ni,nj)
1837  ELSEIF (kj == 3) THEN
1838  mat_loc1(ix,jx) = sij(5,ni,nj)
1839  mat_loc2(ix,jx) = sij(5,ni,nj)
1840  ENDIF
1841  END DO
1842  END DO
1843  END DO
1844  CALL matsetvalues(h_p_phi_mat1, n_ws2, idxn(1:n_ws2), 3*n_ws1, jdxn(1:3*n_ws1), &
1845  mat_loc1(1:n_ws2,1:3*n_ws1), add_values, ierr)
1846  CALL matsetvalues(h_p_phi_mat2, n_ws2, idxn(1:n_ws2), 3*n_ws1, jdxn(1:3*n_ws1), &
1847  mat_loc2(1:n_ws2,1:3*n_ws1), add_values, ierr)
1848 
1849  !Feb 2 2007
1850  mat_loc1 = 0.d0
1851  mat_loc2 = 0.d0
1852  sij = c_sym*sij !SYM
1853  DO ki =1,3
1854  DO ni = 1, n_ws1
1855  i = interface_h_phi%jjs1(ni,ms)
1856  ib = la_h%loc_to_glob(ki,i)
1857  ix = (ki-1)*n_ws1 + ni
1858  idxn(ix) = ib - 1
1859  DO nj = 1, n_ws2
1860  j = interface_h_phi%jjs2(nj,ms)
1861  jb = la_phi%loc_to_glob(1,j)
1862  jx = nj
1863  jdxn(jx) = jb - 1
1864  IF (ki == 1) THEN
1865  mat_loc1(ix,jx) = sij(1,nj,ni)
1866  mat_loc2(ix,jx) = sij(1,nj,ni)
1867  ELSEIF (ki == 2) THEN
1868  mat_loc1(ix,jx) = sij(3,nj,ni)
1869  mat_loc2(ix,jx) = sij(4,nj,ni)
1870  ELSEIF (ki == 3) THEN
1871  mat_loc1(ix,jx) = sij(5,nj,ni)
1872  mat_loc2(ix,jx) = sij(5,nj,ni)
1873  ENDIF
1874  END DO
1875  END DO
1876  END DO
1877  CALL matsetvalues(h_p_phi_mat1, 3*n_ws1, idxn(1:3*n_ws1), n_ws2, jdxn(1:n_ws2), &
1878  mat_loc1(1:3*n_ws1,1:n_ws2), add_values, ierr)
1879  CALL matsetvalues(h_p_phi_mat2, 3*n_ws1, idxn(1:3*n_ws1), n_ws2, jdxn(1:n_ws2), &
1880  mat_loc2(1:3*n_ws1,1:n_ws2), add_values, ierr)
1881  !Feb 2 2007
1882 
1883  sij = 0.d0
1884 
1885  DO ls = 1, l_gs
1886 
1887  !===Compute radius of Gauss point
1888  ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1889 !!$ x = rjs(ls,ms2)*ray/sigma(m1)
1890 ! TEST DEBUG
1891  x = rjs(ls,ms2)*ray/sum(sigma_np(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
1892 ! TEST DEBUG
1893  !terme avec derivee de bi seulement
1894  DO ni = 1, n_ws2
1895  y = x*wws(ni,ls)*mode/ray
1896  DO nj = 1, n_w1
1897  sij(3,ni,nj) = sij(3,ni,nj) + &
1898  y*(dw_cs(2,nj,ls,ms1)*rnorms(2,ls,ms2) + dw_cs(1,nj,ls,ms1)*rnorms(1,ls,ms2))
1899  ENDDO
1900  ENDDO
1901  ENDDO
1902  sij(4,:,:) = -sij(3,:,:)
1903  !TEST
1904  !Sij = 0.d0
1905  !TEST
1906  kj=2
1907  DO ni = 1, n_ws2
1908  i = interface_h_phi%jjs2(ni,ms)
1909  ib = la_phi%loc_to_glob(1,i)
1910  idxn(ni) = ib - 1
1911  DO nj = 1, n_w1
1912  j = h_mesh%jj(nj,m1)
1913  jb = la_h%loc_to_glob(kj,j)
1914  jdxn(nj) = jb - 1
1915  END DO
1916  END DO
1917  CALL matsetvalues(h_p_phi_mat1, n_ws2, idxn(1:n_ws2), n_w1, jdxn(1:n_w1), &
1918  sij(3,1:n_ws2,1:n_w1), add_values, ierr)
1919  CALL matsetvalues(h_p_phi_mat2, n_ws2, idxn(1:n_ws2), n_w1, jdxn(1:n_w1), &
1920  sij(4,1:n_ws2,1:n_w1), add_values, ierr)
1921 
1922  !Feb 2 2007
1923  sij = c_sym*sij !SYM
1924 
1925  mat_loc1 = 0.d0
1926  mat_loc2 = 0.d0
1927  ki=2
1928  DO ni = 1, n_w1
1929  i = h_mesh%jj(ni,m1)
1930  ib = la_h%loc_to_glob(ki,i)
1931  idxn(ni) = ib - 1
1932  DO nj = 1, n_ws2
1933  j = interface_h_phi%jjs2(nj,ms)
1934  jb = la_phi%loc_to_glob(1,j)
1935  jdxn(nj) = jb - 1
1936  mat_loc1(ix,jx) = sij(3,nj,ni)
1937  mat_loc2(ix,jx) = sij(4,nj,ni)
1938  END DO
1939  END DO
1940  CALL matsetvalues(h_p_phi_mat1, n_w1, idxn(1:n_w1), n_ws2, jdxn(1:n_ws2), &
1941  mat_loc1(1:n_w1,1:n_ws2), add_values, ierr)
1942  CALL matsetvalues(h_p_phi_mat2, n_w1, idxn(1:n_w1), n_ws2, jdxn(1:n_ws2), &
1943  mat_loc2(1:n_w1,1:n_ws2), add_values, ierr)
1944  !Feb 2 2007
1945 
1946  sij = 0.d0
1947  DO ls = 1, l_gs
1948 
1949  !===Compute radius of Gauss point
1950  ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1951 !!$ x = rjs(ls,ms2)*ray/sigma(m1)
1952 ! TEST DEBUG
1953  x = rjs(ls,ms2)*ray/sum(sigma_np(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
1954 ! TEST DEBUG
1955  !terme avec derivee de phi et derivee de bi
1956  DO ni = 1, n_w2
1957  y = x*(dw_s(2,ni,ls,ms2)*rnorms(1,ls,ms2) - dw_s(1,ni,ls,ms2)*rnorms(2,ls,ms2))
1958  DO nj = 1, n_w1
1959  sij(1,ni,nj) = sij(1,ni,nj) + y *dw_cs(2,nj,ls,ms1)
1960  sij(5,ni,nj) = sij(5,ni,nj) + (-y)*dw_cs(1,nj,ls,ms1)
1961  ENDDO
1962  ENDDO
1963 
1964  ENDDO
1965 
1966  !TEST
1967  !Sij = 0.d0
1968  !TEST
1969  mat_loc1 = 0.d0
1970  mat_loc2 = 0.d0
1971  DO ni = 1, n_w2
1972  i = phi_mesh%jj(ni,m2)
1973  ib = la_phi%loc_to_glob(1,i)
1974  ix = ni
1975  idxn(ix) = ib - 1
1976  DO nj = 1, n_w1
1977  j = h_mesh%jj(nj,m1)
1978  DO kj=1,3
1979  jb = la_h%loc_to_glob(kj,j)
1980  jx = (kj-1)*n_w1 + nj
1981  jdxn(jx) = jb - 1
1982  IF (kj == 1) THEN
1983  mat_loc1(ix,jx) = sij(1,ni,nj)
1984  mat_loc2(ix,jx) = sij(1,ni,nj)
1985  ELSEIF (kj == 3) THEN
1986  mat_loc1(ix,jx) = sij(5,ni,nj)
1987  mat_loc2(ix,jx) = sij(5,ni,nj)
1988  ENDIF
1989  END DO
1990  END DO
1991  END DO
1992  CALL matsetvalues(h_p_phi_mat1, n_w2, idxn(1:n_w2), 3*n_w1, jdxn(1:3*n_w1), &
1993  mat_loc1(1:n_w2,1:3*n_w1), add_values, ierr)
1994  CALL matsetvalues(h_p_phi_mat2, n_w2, idxn(1:n_w2), 3*n_w1, jdxn(1:3*n_w1), &
1995  mat_loc2(1:n_w2,1:3*n_w1), add_values, ierr)
1996 
1997  !Feb 2 2007
1998  mat_loc1 = 0.d0
1999  mat_loc2 = 0.d0
2000  sij=c_sym*sij !SYM
2001  DO ki=1,3
2002  DO ni = 1, n_w1
2003  i = h_mesh%jj(ni,m1)
2004  ib = la_h%loc_to_glob(ki,i)
2005  ix = (ki-1)*n_w1 + ni
2006  idxn(ix) = ib - 1
2007  DO nj = 1, n_w2
2008  j = phi_mesh%jj(nj,m2)
2009  jb = la_phi%loc_to_glob(1,j)
2010  jx = nj
2011  jdxn(jx) = jb - 1
2012  IF (ki == 1) THEN
2013  mat_loc1(ix,jx) = sij(1,nj,ni)
2014  mat_loc2(ix,jx) = sij(1,nj,ni)
2015  ELSEIF (ki == 3) THEN
2016  mat_loc1(ix,jx) = sij(5,nj,ni)
2017  mat_loc2(ix,jx) = sij(5,nj,ni)
2018  ENDIF
2019  END DO
2020  END DO
2021  END DO
2022  CALL matsetvalues(h_p_phi_mat1, 3*n_w1, idxn(1:3*n_w1), n_w2, jdxn(1:n_w2), &
2023  mat_loc1(1:3*n_w1,1:n_w2), add_values, ierr)
2024  CALL matsetvalues(h_p_phi_mat2, 3*n_w1, idxn(1:3*n_w1), n_w2, jdxn(1:n_w2), &
2025  mat_loc2(1:3*n_w1,1:n_w2), add_values, ierr)
2026 
2027  !JLG, FL, May, 28, 2009
2028  !Ajout du laplacien de phi
2029  sij = 0.d0
2030  DO ls = 1, l_gs
2031  !===Compute radius of Gauss point
2032  ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
2033  muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
2034  x = c_lap*muhl*rjs(ls,ms2)*ray
2035  DO ni = 1, n_ws2
2036  DO nj = 1, n_ws1
2037  sij(1,ni,nj) = sij(1,ni,nj) - x*w_cs(nj,ls)*wws(ni,ls)*rnorms(1,ls,ms2)
2038  sij(5,ni,nj) = sij(5,ni,nj) - x*w_cs(nj,ls)*wws(ni,ls)*rnorms(2,ls,ms2)
2039  ENDDO
2040  END DO
2041  END DO
2042 
2043  mat_loc1 = 0.d0
2044  mat_loc2 = 0.d0
2045  DO ni = 1, n_ws2
2046  i = interface_h_phi%jjs2(ni,ms)
2047  ib = la_phi%loc_to_glob(1,i)
2048  ix = ni
2049  idxn(ix) = ib - 1
2050  DO nj = 1, n_ws1
2051  j = interface_h_phi%jjs1(nj,ms)
2052  jb = la_h%loc_to_glob(1,j)
2053  jx = nj !(1-1)*n_ws1 + nj
2054  jdxn(jx) = jb - 1
2055  mat_loc1(ix,jx) = sij(1,ni,nj)
2056  mat_loc2(ix,jx) = sij(1,ni,nj)
2057 
2058  jb = la_h%loc_to_glob(3,j)
2059  jx = n_ws1 + nj !(3-1)*n_ws1 + nj
2060  jdxn(jx) = jb - 1
2061  mat_loc1(ix,jx) = sij(5,ni,nj)
2062  mat_loc2(ix,jx) = sij(5,ni,nj)
2063  END DO
2064  END DO
2065  CALL matsetvalues(h_p_phi_mat1, n_ws2, idxn(1:n_ws2), 2*n_ws1, jdxn(1:2*n_ws1), &
2066  mat_loc1(1:n_ws2,1:2*n_ws1), add_values, ierr)
2067  CALL matsetvalues(h_p_phi_mat2, n_ws2, idxn(1:n_ws2), 2*n_ws1, jdxn(1:2*n_ws1), &
2068  mat_loc2(1:n_ws2,1:2*n_ws1), add_values, ierr)
2069  !JLG, FL, May, 28, 2009
2070 
2071  !Feb 2 2007
2072  !==================
2073 
2074  !(use .true. for convergence tests)
2075  !June 6 2008, I put back (.true.) always.
2076  !Works much better when mu is discontinuous.
2077  !Mars 22 2007
2078 
2079  IF (stab(2) > 1.d-12) THEN
2080  !IF (.FALSE.) THEN
2081  !Mars 22 2007
2082  !Enforcing weak continuity on the normal components
2083  hsij = 0.d0
2084  sij = 0.d0
2085  phisij = 0.d0
2086 
2087 
2088  ms2 = interface_h_phi%mesh2(ms)
2089  hm1 = sum(rjs(:,ms2))**(2*alpha-1)
2090 
2091  DO ls = 1, l_gs
2092 
2093  !Feb 8 2007, muhl
2094  muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
2095  !Feb 8 2007, muhl
2096  ray = 0.d0
2097  DO ni = 1, n_ws2; i = phi_mesh%jjs(ni,ms2)
2098  ray = ray + phi_mesh%rr(1,i)* phi_mesh%gauss%wws(ni,ls)
2099  END DO
2100 
2101 
2102  !ray = ray*hm1*rjs(ls,ms2)
2103  !June 8, 2008, Normalization, JLG, FL, May, 28, 2009
2104  ray = stab_div*ray*hm1*rjs(ls,ms2)
2105  !ray = stab_div*ray*hm1*rjs(ls,ms2)/muhl
2106  !ray = stab_div*ray*hm1*rjs(ls,ms2)/muhl**2
2107  !June 8, 2008, Normalization, JLG, FL, May, 28, 2009
2108  DO ni = 1, n_ws1
2109  DO nj = 1, n_ws1
2110  x = muhl**2*w_cs(ni,ls)*w_cs(nj,ls)*ray
2111  hsij(1,ni,nj) = hsij(1,ni,nj) + x*rnorms(1,ls,ms2)**2
2112  hsij(4,ni,nj) = hsij(4,ni,nj) + x*rnorms(1,ls,ms2)*rnorms(2,ls,ms2)
2113  hsij(6,ni,nj) = hsij(6,ni,nj) + x*rnorms(2,ls,ms2)**2
2114  END DO
2115 
2116  DO nj = 1, n_w2
2117  x = muhl*mu_phi*w_cs(ni,ls)*(dw_s(1,nj,ls,ms2)*rnorms(1,ls,ms2) +&
2118  dw_s(2,nj,ls,ms2)*rnorms(2,ls,ms2))*ray
2119  sij(1,ni,nj) = sij(1,ni,nj) - x*rnorms(1,ls,ms2)
2120  sij(5,ni,nj) = sij(5,ni,nj) - x*rnorms(2,ls,ms2)
2121  ENDDO
2122  ENDDO
2123 
2124  DO ni = 1, n_w2
2125  DO nj = 1, n_w2
2126  x = mu_phi**2*(dw_s(1,ni,ls,ms2)*rnorms(1,ls,ms2) + dw_s(2,ni,ls,ms2)*rnorms(2,ls,ms2))* &
2127  (dw_s(1,nj,ls,ms2)*rnorms(1,ls,ms2) + dw_s(2,nj,ls,ms2)*rnorms(2,ls,ms2))*ray
2128  phisij(ni,nj) = phisij(ni,nj) + x
2129  ENDDO
2130  ENDDO
2131 
2132  END DO
2133  sij(2,:,:) = sij(1,:,:)
2134  sij(6,:,:) = sij(5,:,:)
2135 
2136 
2137  mat_loc1 = 0.d0
2138  mat_loc2 = 0.d0
2139  DO ni = 1, n_ws1
2140  i = h_mesh%jjs(ni,ms1)
2141  DO ki= 1, 3, 2
2142  ib = la_h%loc_to_glob(ki,i)
2143  ix = (ki/2)*n_ws1 + ni
2144  idxn(ix) = ib - 1
2145  DO nj = 1, n_ws1
2146  j = h_mesh%jjs(nj,ms1)
2147  DO kj = 1, 3, 2
2148  jb = la_h%loc_to_glob(kj,j)
2149  jx = (kj/2)*n_ws1 + nj
2150  jdxn(jx) = jb - 1
2151  IF (ki*kj==1) THEN
2152  mat_loc1(ix,jx) = hsij(1,ni,nj)
2153  mat_loc2(ix,jx) = hsij(1,ni,nj)
2154  ELSE IF (ki*kj==9) THEN
2155  mat_loc1(ix,jx) = hsij(6,ni,nj)
2156  mat_loc2(ix,jx) = hsij(6,ni,nj)
2157  ELSE IF (ki*kj==3) THEN
2158  mat_loc1(ix,jx) = hsij(4,ni,nj)
2159  mat_loc2(ix,jx) = hsij(4,ni,nj)
2160  END IF
2161  END DO
2162  END DO
2163 
2164  DO nj = 1, n_w2
2165  j = phi_mesh%jj(nj,m2)
2166  jb = la_phi%loc_to_glob(1,j)
2167  jx = 2*n_ws1 + nj
2168  jdxn(jx) = jb - 1
2169  mat_loc1(ix,jx) = sij(2*ki-1,ni,nj)
2170  mat_loc2(ix,jx) = sij(2*ki-1,ni,nj)
2171  END DO
2172  ENDDO
2173  ENDDO
2174  CALL matsetvalues(h_p_phi_mat1, 2*n_ws1, idxn(1:2*n_ws1), 2*n_ws1+n_w2, jdxn(1:2*n_ws1+n_w2), &
2175  mat_loc1(1:2*n_ws1,1:2*n_ws1+n_w2), add_values, ierr)
2176  CALL matsetvalues(h_p_phi_mat2, 2*n_ws1, idxn(1:2*n_ws1), 2*n_ws1+n_w2, jdxn(1:2*n_ws1+n_w2), &
2177  mat_loc2(1:2*n_ws1,1:2*n_ws1+n_w2), add_values, ierr)
2178 
2179  mat_loc1 = 0.d0
2180  mat_loc2 = 0.d0
2181  DO ni = 1, n_w2
2182  i = phi_mesh%jj(ni,m2)
2183  ib = la_phi%loc_to_glob(1,i)
2184  ix = ni
2185  idxn(ix) = ib -1
2186  DO nj = 1, n_ws1
2187  j = h_mesh%jjs(nj,ms1)
2188  DO kj = 1, 3, 2
2189  jb = la_h%loc_to_glob(kj,j)
2190  jx = (kj/2)*n_ws1 + nj
2191  jdxn(jx) = jb - 1
2192  mat_loc1(ix,jx) = sij(2*kj-1,nj,ni)
2193  mat_loc2(ix,jx) = sij(2*kj-1,nj,ni)
2194  END DO
2195  END DO
2196 
2197  DO nj = 1, n_w2
2198  j = phi_mesh%jj(nj,m2)
2199  jb = la_phi%loc_to_glob(1,j)
2200  jx = 2*n_ws1 + nj
2201  jdxn(jx) = jb - 1
2202  mat_loc1(ix,jx) = phisij(ni,nj)
2203  mat_loc2(ix,jx) = phisij(ni,nj)
2204  END DO
2205  END DO
2206  CALL matsetvalues(h_p_phi_mat1, n_w2, idxn(1:n_w2), 2*n_ws1+n_w2, jdxn(1:2*n_ws1+n_w2), &
2207  mat_loc1(1:n_w2,1:2*n_ws1+n_w2), add_values, ierr)
2208  CALL matsetvalues(h_p_phi_mat2, n_w2, idxn(1:n_w2), 2*n_ws1+n_w2, jdxn(1:2*n_ws1+n_w2), &
2209  mat_loc2(1:n_w2,1:2*n_ws1+n_w2), add_values, ierr)
2210  END IF
2211  !FIN TEST
2212 
2213  ENDDO
2214 
2215 
2216  !=========================================================
2217  !--- Artificial boundary condition: d(phi)/dR + (1/R)*phi = 0
2218  !=========================================================
2219 
2220  IF (.NOT.present(index_fourier) .OR. .NOT.present(r_fourier)) RETURN
2221  IF (r_fourier.GT.0.d0) THEN
2222  !WRITE(*,*) ' Assembling the Fourier condition'
2223  DO ms = 1, phi_mesh%mes
2224  IF (phi_mesh%sides(ms) /= index_fourier) cycle ! Not on the artificial boundary
2225 
2226  phisij = 0.d0
2227 
2228  DO ls = 1, phi_mesh%gauss%l_Gs
2229 
2230  !===Compute radius of Gauss point
2231  ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms))* phi_mesh%gauss%wws(:,ls))
2232 
2233  x = c_mu_phi*rjs(ls,ms)*ray/r_fourier
2234 
2235  DO ni=1, phi_mesh%gauss%n_ws
2236  DO nj=1, phi_mesh%gauss%n_ws
2237  phisij(ni,nj) = phisij(ni,nj) + x*wws(ni,ls)*wws(nj,ls)
2238  ENDDO
2239  ENDDO
2240 
2241  ENDDO
2242 
2243 
2244  DO ni = 1, phi_mesh%gauss%n_ws
2245  i = phi_mesh%jjs(ni,ms)
2246  ib = la_phi%loc_to_glob(1,i)
2247  idxn(ni) = ib - 1
2248  DO nj = 1, phi_mesh%gauss%n_ws
2249  j = phi_mesh%jjs(nj,ms)
2250  jb = la_phi%loc_to_glob(1,j)
2251  jdxn(nj) = jb - 1
2252  END DO
2253  END DO
2254  CALL matsetvalues(h_p_phi_mat1, n_ws2, idxn(1:n_ws2), n_ws2, jdxn(1:n_ws2), &
2255  phisij(1:n_ws2,1:n_ws2), add_values, ierr)
2256  CALL matsetvalues(h_p_phi_mat2, n_ws2, idxn(1:n_ws2), n_ws2, jdxn(1:n_ws2), &
2257  phisij(1:n_ws2,1:n_ws2), add_values, ierr)
2258  END DO
2259  END IF
2260 
2261  CALL matassemblybegin(h_p_phi_mat1,mat_final_assembly,ierr)
2262  CALL matassemblyend(h_p_phi_mat1,mat_final_assembly,ierr)
2263  CALL matassemblybegin(h_p_phi_mat2,mat_final_assembly,ierr)
2264  CALL matassemblyend(h_p_phi_mat2,mat_final_assembly,ierr)
2265 
2266 !!$ DEALLOCATE(mat_loc1, mat_loc2, idxn, jdxn)
2267 
2268  END SUBROUTINE mat_h_p_phi_maxwell
2269 
2270  SUBROUTINE mat_dirichlet_maxwell(H_mesh, Dirichlet_bdy_H_sides, &
2271  mode, stab, la_h, h_p_phi_mat1, h_p_phi_mat2, sigma_np)
2272  USE def_type_mesh
2273  USE dir_nodes
2274  USE gauss_points
2275  USE boundary
2276  USE my_util
2277  IMPLICIT NONE
2278  TYPE(mesh_type), INTENT(IN) :: h_mesh
2279  INTEGER, DIMENSION(:), INTENT(IN) :: dirichlet_bdy_h_sides
2280  INTEGER, INTENT(IN) :: mode
2281  REAL(KIND=8), DIMENSION(3), INTENT(IN) :: stab
2282  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: sigma_np
2283 
2284  INTEGER :: ms, ls, ni, nj, i, j, &
2285  n_ws1, n_w1, m1, ki, kj, ib, jb
2286  REAL(KIND=8) :: x, y, hm1
2287  REAL(KIND=8) :: ray, error, stab_colle_h_mu
2288  REAL(KIND=8), DIMENSION(9,H_mesh%gauss%n_w,H_mesh%gauss%n_w) :: hsij
2289  ! MATRICES POUR LES TERMES DE BORDS Hsij et Phisij
2290  !=================================================
2291  ! (--------------------------------------------------------------------)
2292  ! ( Hsij(1) | Hsij(2) | Hsij(4) || Sij(1) )
2293  ! ( Hsij(1) | Hsij(3) | Hsij(4) || Sij(2) )
2294  ! (--------------------------------------------------------------------)
2295  ! ( | Hsij(5) | || Sij(3) )
2296  ! ( | Hsij(5) | || Sij(4) )
2297  ! (--------------------------------------------------------------------)
2298  ! ( Hsij(7) | Hsij(9) | Hsij(6) || Sij(5) )
2299  ! ( Hsij(7) | Hsij(8) | Hsij(6) || Sij(6) )
2300  ! (====================================================================)
2301  ! ( Sij'(1) | Sij'(3) | Sij'(5) || Phisij )
2302  ! ( Sij'(2) | Sij'(4) | Sij'(6) || Phisij )
2303  ! (------------------------------------------------------------------- )
2304  !
2305  ! L'autre partie des termes croises est la symetrique de la premiere
2306  ! juste apres le calcsrhs_maul du terme de bord dissymetrique
2307  !June 8 2008
2308  REAL(KIND=8) :: c_sym=.0d0 ! Symmetrization of the bilinear form
2309  !June 8 2008
2310 !!$ FL + CN 22/03 2013
2311 !!$ REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: mat_loc1, mat_loc2
2312 !!$ INTEGER , DIMENSION(:), ALLOCATABLE :: idxn, jdxn
2313  REAL(KIND=8), DIMENSION(3*H_mesh%gauss%n_w,3*H_mesh%gauss%n_w) :: mat_loc1, mat_loc2
2314  INTEGER , DIMENSION(3*H_mesh%gauss%n_w) :: idxn, jdxn
2315 !!$ FL + CN 22/03 2013
2316  TYPE(petsc_csr_la) :: la_h
2317  INTEGER :: ix, jx
2318  INTEGER :: count
2319 #include "petsc/finclude/petsc.h"
2320  petscerrorcode :: ierr
2321  mat :: h_p_phi_mat1, h_p_phi_mat2
2322 
2323  !June 2009, JLG, CN, Normalization
2324  stab_colle_h_mu = stab(3)
2325  !Jan 2010, JLG, CN, Normalization,
2326 
2327  !*********************************************************************************
2328  !--------------------TERMS ON DIRICHLET BOUNDARY-----------------------------
2329  !**********************************************************************************
2330  CALL gauss(h_mesh)
2331  n_ws1 = h_mesh%gauss%n_ws
2332  n_w1 = h_mesh%gauss%n_w
2333 
2334 !!$ ALLOCATE(mat_loc1(3*n_w1,3*n_w1))
2335 !!$ ALLOCATE(mat_loc2(3*n_w1,3*n_w1))
2336 !!$ ALLOCATE(idxn(3*n_w1))
2337 !!$ ALLOCATE(jdxn(3*n_w1))
2338 
2339  error = 0
2340  DO count = 1, SIZE(dirichlet_bdy_h_sides)
2341  ms = dirichlet_bdy_h_sides(count)
2342  hm1 = stab_colle_h_mu/sum(h_mesh%gauss%rjs(:,ms))
2343  m1 = h_mesh%neighs(ms)
2344  !====================================================================================
2345  !------------------------------------TERMES SUR LE BLOC H----------------------------
2346  !====================================================================================
2347 
2348  !-------------------------------hm1 (bi x ni) . (bj x nj)----------------------------
2349  !====================================================================================
2350 
2351  hsij = 0.d0
2352  DO ls = 1, h_mesh%gauss%l_Gs
2353  !===Compute radius of Gauss point
2354  ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
2355  x = hm1*h_mesh%gauss%rjs(ls,ms)*ray
2356 
2357  DO ni = 1, h_mesh%gauss%n_ws
2358  DO nj = 1, h_mesh%gauss%n_ws
2359  y = x * h_mesh%gauss%wws(ni,ls)*h_mesh%gauss%wws(nj,ls)
2360  hsij(1,ni,nj) = hsij(1,ni,nj) + y*(h_mesh%gauss%rnorms(2,ls,ms)**2)
2361  hsij(4,ni,nj) = hsij(4,ni,nj) - y*h_mesh%gauss%rnorms(1,ls,ms)*h_mesh%gauss%rnorms(2,ls,ms)
2362  hsij(5,ni,nj) = hsij(5,ni,nj) + y
2363  hsij(6,ni,nj) = hsij(6,ni,nj) + y*(h_mesh%gauss%rnorms(1,ls,ms)**2)
2364  ENDDO
2365  ENDDO
2366 
2367  ENDDO
2368 
2369 
2370  !TEST
2371  !Hsij = 0.d0
2372  !Hsij = Hsij / hm1
2373  !TEST
2374  mat_loc1 = 0.d0
2375  mat_loc2 = 0.d0
2376  DO ki= 1, 3
2377  DO ni = 1, n_ws1
2378  i = h_mesh%jjs(ni,ms)
2379  ib = la_h%loc_to_glob(ki,i)
2380  ix = ni + (ki-1)*n_ws1
2381  idxn(ix) = ib - 1
2382  DO kj = 1, 3
2383  DO nj = 1, n_ws1
2384  j = h_mesh%jjs(nj,ms)
2385  jb = la_h%loc_to_glob(kj,j)
2386  jx = nj + (kj-1)*n_ws1
2387  jdxn(jx) = jb - 1
2388  IF ((ki == 1) .AND. (kj == 1)) THEN
2389  mat_loc1(ix,jx) = hsij(1,ni,nj)
2390  mat_loc2(ix,jx) = hsij(1,ni,nj)
2391  ELSEIF ((ki == 1) .AND. (kj == 3)) THEN
2392  mat_loc1(ix,jx) = hsij(4,ni,nj)
2393  mat_loc2(ix,jx) = hsij(4,ni,nj)
2394  ELSEIF ((ki == 3) .AND. (kj == 1)) THEN
2395  mat_loc1(ix,jx) = hsij(4,nj,ni)
2396  mat_loc2(ix,jx) = hsij(4,nj,ni)
2397  ELSEIF ((ki == 2) .AND. (kj == 2)) THEN
2398  mat_loc1(ix,jx) = hsij(5,ni,nj)
2399  mat_loc2(ix,jx) = hsij(5,ni,nj)
2400  ELSEIF ((ki == 3) .AND. (kj == 3)) THEN
2401  mat_loc1(ix,jx) = hsij(6,ni,nj)
2402  mat_loc2(ix,jx) = hsij(6,ni,nj)
2403  ENDIF
2404  END DO
2405  END DO
2406  END DO
2407  END DO
2408 
2409  CALL matsetvalues(h_p_phi_mat1, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
2410  mat_loc1(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
2411  CALL matsetvalues(h_p_phi_mat2, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
2412  mat_loc2(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
2413 
2414  !====================================================================================
2415  !------------------------(1/sigma) (Rot bj) . (bi x ni)------------------------------
2416  !====================================================================================
2417 
2418  !JLG+FL: Jan 18 2013
2419  !There was a bug on the sign of the normal
2420  !The sign before rnorms has been changed everywhere in this loop.
2421  hsij = 0.d0
2422 
2423  DO ls = 1, h_mesh%gauss%l_Gs
2424  !===Compute radius of Gauss point
2425  ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms))* h_mesh%gauss%wws(:,ls))
2426  x = h_mesh%gauss%rjs(ls,ms)*ray/sum(sigma_np(h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
2427 
2428  !terme sans derivees
2429  DO ni = 1,n_ws1
2430  DO nj = 1, n_ws1
2431  y = x*h_mesh%gauss%wws(ni,ls)*h_mesh%gauss%wws(nj,ls)
2432  hsij(2,ni,nj) = hsij(2,ni,nj) + y * (-mode/ray)*(rnorms(1,ls,ms))
2433  hsij(3,ni,nj) = hsij(3,ni,nj) + y * mode/ray *(rnorms(1,ls,ms))
2434  hsij(5,ni,nj) = hsij(5,ni,nj) + y * (-1/ray) *(rnorms(1,ls,ms))
2435  hsij(8,ni,nj) = hsij(8,ni,nj) + y * (-mode/ray)*(rnorms(2,ls,ms))
2436  hsij(9,ni,nj) = hsij(9,ni,nj) + y * mode/ray *(rnorms(2,ls,ms))
2437  ENDDO
2438  ENDDO
2439 
2440  ENDDO
2441 
2442  !TEST
2443  !Hsij = 0.d0
2444  !TEST
2445 
2446  mat_loc1 = 0.d0
2447  mat_loc2 = 0.d0
2448  DO ki= 1, 3
2449  DO ni = 1, n_ws1
2450  i = h_mesh%jjs(ni,ms)
2451  ib = la_h%loc_to_glob(ki,i)
2452  ix = ni + (ki-1)*n_ws1
2453  idxn(ix) = ib - 1
2454  DO kj = 1, 3
2455  DO nj = 1, n_ws1
2456  j = h_mesh%jjs(nj,ms)
2457  jb = la_h%loc_to_glob(kj,j)
2458  jx = nj + (kj-1)*n_ws1
2459  jdxn(jx) = jb - 1
2460  IF ( (ki == 2) .AND. (kj == 1)) THEN
2461  mat_loc1(ix,jx) = hsij(2,ni,nj)
2462  mat_loc2(ix,jx) = hsij(3,ni,nj)
2463  ELSEIF ((ki == 2) .AND. (kj == 2)) THEN
2464  mat_loc1(ix,jx) = hsij(5,ni,nj)
2465  mat_loc2(ix,jx) = hsij(5,ni,nj)
2466  ELSEIF ( (ki == 2) .AND. (kj == 3)) THEN
2467  mat_loc1(ix,jx) = hsij(8,ni,nj)
2468  mat_loc2(ix,jx) = hsij(9,ni,nj)
2469  ENDIF
2470  END DO
2471  END DO
2472  END DO
2473  END DO
2474 
2475  CALL matsetvalues(h_p_phi_mat1, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
2476  mat_loc1(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
2477  CALL matsetvalues(h_p_phi_mat2, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
2478  mat_loc2(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
2479 
2480  !Feb 2 2007
2481  hsij=c_sym*hsij !SYM
2482  mat_loc1 = 0.d0
2483  mat_loc2 = 0.d0
2484  DO ki= 1, 3
2485  DO ni = 1, n_ws1
2486  i = h_mesh%jjs(ni,ms)
2487  ib = la_h%loc_to_glob(ki,i)
2488  ix = ni + (ki-1)*n_ws1
2489  idxn(ix) = ib - 1
2490  DO kj = 1, 3
2491  DO nj = 1, n_ws1
2492  j = h_mesh%jjs(nj,ms)
2493  jb = la_h%loc_to_glob(kj,j)
2494  jx = nj + (kj-1)*n_ws1
2495  jdxn(jx) = jb - 1
2496  IF ( (kj == 2) .AND. (ki == 1)) THEN
2497  mat_loc1(ix,jx) = hsij(2,nj,ni)
2498  mat_loc2(ix,jx) = hsij(3,nj,ni)
2499  ELSEIF ((kj == 2) .AND. (ki == 2)) THEN
2500  mat_loc1(ix,jx) = hsij(5,nj,ni)
2501  mat_loc2(ix,jx) = hsij(5,nj,ni)
2502  ELSEIF ( (kj == 2) .AND. (ki == 3)) THEN
2503  mat_loc1(ix,jx) = hsij(8,nj,ni)
2504  mat_loc2(ix,jx) = hsij(9,nj,ni)
2505  ENDIF
2506  END DO
2507  END DO
2508  END DO
2509  END DO
2510  CALL matsetvalues(h_p_phi_mat1, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
2511  mat_loc1(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
2512  CALL matsetvalues(h_p_phi_mat2, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
2513  mat_loc2(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
2514  !feb 2 2007
2515 
2516 
2517  hsij = 0.d0
2518 
2519  DO ls = 1, h_mesh%gauss%l_Gs
2520  !===Compute radius of Gauss point
2521  ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms))* h_mesh%gauss%wws(:,ls))
2522  x = h_mesh%gauss%rjs(ls,ms)*ray /sum(sigma_np(h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
2523 
2524  !termes avec derivees
2525  DO ni = 1,n_ws1
2526  y = x*h_mesh%gauss%wws(ni,ls)
2527  DO nj = 1, n_w1
2528  hsij(1,ni,nj) = hsij(1,ni,nj) + y*(-h_mesh%gauss%dw_s(2,nj,ls,ms))*(rnorms(2,ls,ms))
2529  hsij(4,ni,nj) = hsij(4,ni,nj) + y* h_mesh%gauss%dw_s(1,nj,ls,ms) *(rnorms(2,ls,ms))
2530  hsij(5,ni,nj) = hsij(5,ni,nj) + &
2531  y*(-h_mesh%gauss%dw_s(2,nj,ls,ms)*(rnorms(2,ls,ms))-h_mesh%gauss%dw_s(1,nj,ls,ms)*(rnorms(1,ls,ms)))
2532  hsij(6,ni,nj) = hsij(6,ni,nj) + y*(-h_mesh%gauss%dw_s(1,nj,ls,ms))*(rnorms(1,ls,ms))
2533  hsij(7,ni,nj) = hsij(7,ni,nj) + y* h_mesh%gauss%dw_s(2,nj,ls,ms) *(rnorms(1,ls,ms))
2534  ENDDO
2535  ENDDO
2536 
2537  ENDDO
2538 
2539  !TEST
2540  !Hsij = 0.d0
2541  !TEST
2542 
2543 
2544  mat_loc1 = 0.d0
2545  mat_loc2 = 0.d0
2546  DO ki= 1, 3
2547  DO ni = 1, n_ws1
2548  i = h_mesh%jjs(ni,ms)
2549  ib = la_h%loc_to_glob(ki,i)
2550  ix = ni + (ki-1)*n_ws1
2551  idxn(ix) = ib - 1
2552  DO kj = 1, 3
2553  DO nj = 1, n_w1
2554  j = h_mesh%jj(nj,m1)
2555  jb = la_h%loc_to_glob(kj,j)
2556  jx = nj + (kj-1)*n_w1
2557  jdxn(jx) = jb - 1
2558  IF ((ki == 1) .AND. (kj == 1)) THEN
2559  mat_loc1(ix,jx) = hsij(1,ni,nj)
2560  mat_loc2(ix,jx) = hsij(1,ni,nj)
2561  ELSEIF ((ki == 1) .AND. (kj == 3)) THEN
2562  mat_loc1(ix,jx) = hsij(4,ni,nj)
2563  mat_loc2(ix,jx) = hsij(4,ni,nj)
2564  ELSEIF ((ki == 2) .AND. (kj == 2)) THEN
2565  mat_loc1(ix,jx) = hsij(5,ni,nj)
2566  mat_loc2(ix,jx) = hsij(5,ni,nj)
2567  ELSEIF ((ki == 3) .AND. (kj == 3)) THEN
2568  mat_loc1(ix,jx) = hsij(6,ni,nj)
2569  mat_loc2(ix,jx) = hsij(6,ni,nj)
2570  ELSEIF ((ki == 3) .AND. (kj == 1)) THEN
2571  mat_loc1(ix,jx) = hsij(7,ni,nj)
2572  mat_loc2(ix,jx) = hsij(7,ni,nj)
2573  ENDIF
2574  END DO
2575  END DO
2576  END DO
2577  END DO
2578 
2579  CALL matsetvalues(h_p_phi_mat1, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_w1, jdxn(1:3*n_w1), &
2580  mat_loc1(1:3*n_ws1,1:3*n_w1), add_values, ierr)
2581  CALL matsetvalues(h_p_phi_mat2, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_w1, jdxn(1:3*n_w1), &
2582  mat_loc2(1:3*n_ws1,1:3*n_w1), add_values, ierr)
2583 
2584  !Feb 2 2007
2585  hsij=c_sym*hsij !SYM
2586  mat_loc1 = 0.d0
2587  mat_loc2 = 0.d0
2588  DO ki = 1, 3
2589  DO ni = 1, n_w1
2590  i = h_mesh%jj(ni,m1)
2591  ib = la_h%loc_to_glob(ki,i)
2592  ix = ni + (ki-1)*n_w1
2593  idxn(ix) = ib - 1
2594  DO kj= 1, 3
2595  DO nj = 1, n_ws1
2596  j = h_mesh%jjs(nj,ms)
2597  jb = la_h%loc_to_glob(kj,j)
2598  jx = nj + (kj-1)*n_ws1
2599  jdxn(jx) = jb - 1
2600  IF ((kj == 1) .AND. (ki == 1)) THEN
2601  mat_loc1(ix,jx) = hsij(1,nj,ni)
2602  mat_loc2(ix,jx) = hsij(1,nj,ni)
2603  ELSEIF ((kj == 1) .AND. (ki == 3)) THEN
2604  mat_loc1(ix,jx) = hsij(4,nj,ni)
2605  mat_loc2(ix,jx) = hsij(4,nj,ni)
2606  ELSEIF ((kj == 2) .AND. (ki == 2)) THEN
2607  mat_loc1(ix,jx) = hsij(5,nj,ni)
2608  mat_loc2(ix,jx) = hsij(5,nj,ni)
2609  ELSEIF ((kj == 3) .AND. (ki == 3)) THEN
2610  mat_loc1(ix,jx) = hsij(6,nj,ni)
2611  mat_loc2(ix,jx) = hsij(6,nj,ni)
2612  ELSEIF ((kj == 3) .AND. (ki == 1)) THEN
2613  mat_loc1(ix,jx) = hsij(7,nj,ni)
2614  mat_loc2(ix,jx) = hsij(7,nj,ni)
2615  ENDIF
2616  END DO
2617  END DO
2618  END DO
2619  END DO
2620 
2621  CALL matsetvalues(h_p_phi_mat1, 3*n_w1 , idxn(1:3*n_w1), 3*n_ws1, jdxn(1:3*n_ws1), &
2622  mat_loc1(1:3*n_w1,1:3*n_ws1), add_values, ierr)
2623  CALL matsetvalues(h_p_phi_mat2, 3*n_w1 , idxn(1:3*n_w1), 3*n_ws1, jdxn(1:3*n_ws1), &
2624  mat_loc2(1:3*n_w1,1:3*n_ws1), add_values, ierr)
2625 
2626  ENDDO
2627 
2628  CALL matassemblybegin(h_p_phi_mat1,mat_final_assembly,ierr)
2629  CALL matassemblyend(h_p_phi_mat1,mat_final_assembly,ierr)
2630  CALL matassemblybegin(h_p_phi_mat2,mat_final_assembly,ierr)
2631  CALL matassemblyend(h_p_phi_mat2,mat_final_assembly,ierr)
2632 
2633 !!$ IF (ALLOCATED(mat_loc1)) DEALLOCATE(mat_loc1)
2634 !!$ IF (ALLOCATED(mat_loc2)) DEALLOCATE(mat_loc2)
2635 !!$ IF (ALLOCATED(idxn)) DEALLOCATE(idxn)
2636 !!$ IF (ALLOCATED(jdxn)) DEALLOCATE(jdxn)
2637 
2638 
2639  END SUBROUTINE mat_dirichlet_maxwell
2640 
2641  SUBROUTINE courant_int_by_parts(H_mesh,phi_mesh,interface_H_phi,sigma,mu_phi,mu_H_field,time,mode,&
2642  rhs_h,nl, la_h, la_phi, vb_1, vb_2, h_ext, sigma_curl, j_over_sigma)
2643  !forcage faisant intervenir J, volumique et interface_H_phi
2644  !pour le probleme en entier
2645 
2646  USE def_type_mesh
2647  USE gauss_points
2648  USE boundary
2649  USE my_util
2650  USE input_data
2651  IMPLICIT NONE
2652  TYPE(mesh_type), INTENT(IN) :: h_mesh, phi_mesh
2653  TYPE(interface_type), INTENT(IN) :: interface_h_phi
2654  REAL(KIND=8), INTENT(IN) :: mu_phi, time
2655  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: sigma
2656  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: mu_h_field
2657  INTEGER, INTENT(IN) :: mode
2658  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: nl
2659  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: h_ext
2660  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rhs_h
2661  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: j_over_sigma !Used only if sigma variable in fluid
2662  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: sigma_curl !Used only if sigma variable in fluid
2663  REAL(KIND=8), DIMENSION(6) :: j_over_sigma_gauss
2664  INTEGER :: index
2665  REAL(KIND=8), DIMENSION(H_mesh%np,6) :: src_h
2666  REAL(KIND=8), DIMENSION(phi_mesh%np,2) :: src_phi
2667 ! CN possible faute POINTER
2668  !REAL(KIND=8), DIMENSION(:,:), POINTER :: src_H, src_phi
2669  !REAL(KIND=8), DIMENSION(:,:), POINTER :: nl
2670  !REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: src_H, src_phi
2671 !!$ REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: nl
2672  REAL(KIND=8), DIMENSION(phi_mesh%gauss%n_ws,phi_mesh%gauss%l_Gs) :: w_cs
2673  REAL(KIND=8), DIMENSION(2) :: gaussp
2674  REAL(KIND=8) :: ray
2675  INTEGER :: m, l, i, ni, k, ms, ls, n_ws1, n_ws2, ms1, ms2, h_bloc_size, n_w2, m1
2676  INTEGER :: mesh_id1
2677  REAL(KIND=8), DIMENSION(6) :: jsolh_anal, rhs_hl
2678  REAL(KIND=8), DIMENSION(2,H_mesh%gauss%n_w) :: dwh
2679  !REAL(KIND=8), DIMENSION(2,phi_mesh%gauss%n_w) :: dwphi
2680  !REAL(KIND=8), DIMENSION(2,phi_mesh%gauss%n_w) :: src_phil
2681  REAL(KIND=8) :: ray_rjl, muhl
2682  !REAL(KIND=8) :: moderay2
2683  !REAL(KIND=8) :: tps, dummy
2684 !!$ FL + CN, 22/03/2013
2685 !!$ INTEGER, DIMENSION(:), ALLOCATABLE :: idxn
2686  INTEGER, DIMENSION(H_mesh%np) :: idxn_h
2687  INTEGER, DIMENSION(phi_mesh%np) :: idxn_phi
2688 !!$ FL + CN, 22/03/2013
2689  TYPE(petsc_csr_la) :: la_h, la_phi
2690  REAL(KIND=8), DIMENSION(6) :: h_ext_l
2691 #include "petsc/finclude/petsc.h"
2692  petscerrorcode :: ierr
2693  vec :: vb_1, vb_2
2694 
2695  !ALLOCATE(src_H(H_mesh%np,6), src_phi(phi_mesh%np,2))
2696 
2697  CALL veczeroentries(vb_1, ierr)
2698  CALL veczeroentries(vb_2, ierr)
2699 
2700  !forcage volumique
2701  !attention on comprime le calcul sur les points de Gauss et integration !!
2702  !j/sigma *(Rot(b))
2703 
2704  !tps = user_time(dummy)
2705  src_h = 0.d0
2706  src_phi = 0.d0
2707  index = 0
2708 
2709  DO m = 1, h_mesh%me
2710  mesh_id1 = h_mesh%i_d(m)
2711  DO l = 1, h_mesh%gauss%l_G
2712  index = index + 1
2713  !Feb 8 2007, muhl
2714  muhl=sum(mu_h_field(h_mesh%jj(:,m))*h_mesh%gauss%ww(:,l))
2715  !Feb 8 2007, muhl
2716  dwh = h_mesh%gauss%dw(:,:,l,m)
2717  !===Compute radius of Gauss point
2718  DO k=1, 6
2719  h_ext_l(k) = sum(h_ext(h_mesh%jj(:,m),k)*h_mesh%gauss%ww(:,l))
2720  END DO
2721 
2722  jsolh_anal = 0.d0
2723  rhs_hl = 0.d0
2724  gaussp = 0.d0
2725  j_over_sigma_gauss = 0.d0
2726  DO ni = 1, h_mesh%gauss%n_w; i = h_mesh%jj(ni,m)
2727  gaussp = gaussp + h_mesh%rr(:,i)*h_mesh%gauss%ww(ni,l)
2728  jsolh_anal(:) = jsolh_anal(:) + muhl*nl(i,:)*h_mesh%gauss%ww(ni,l)
2729  rhs_hl(:) = rhs_hl(:) + rhs_h(i,:)*h_mesh%gauss%ww(ni,l)
2730  j_over_sigma_gauss(:) = j_over_sigma_gauss(:) + j_over_sigma(i,:)*h_mesh%gauss%ww(ni,l)
2731  ENDDO
2732  ray = gaussp(1)
2733  ray_rjl = h_mesh%gauss%rj(l,m)*ray
2734 
2735  IF (inputs%if_level_set.AND.inputs%variation_sigma_fluid) THEN
2736  DO k = 1, 6
2737  jsolh_anal(k) = jsolh_anal(k) + j_over_sigma_gauss(k) + sigma_curl(index,k)
2738  END DO
2739  ELSE
2740  DO k = 1, 6
2741  !JsolH_anal(k) = muhl*JsolH_anal(k) + & ! BUG Jan 2010, JLG, CN, FL
2742  jsolh_anal(k) = jsolh_anal(k) + &
2743  jexact_gauss(k, gaussp, mode, mu_phi, sigma(m), muhl, time, mesh_id1, h_ext_l)/sigma(m)
2744  END DO
2745  END IF
2746 
2747  DO ni = 1,h_mesh%gauss%n_w
2748 
2749  i = h_mesh%jj(ni,m)
2750 
2751  !--------Composante r------
2752  src_h(i,1) = src_h(i,1)+ ray_rjl &
2753  *(jsolh_anal(3)*dwh(2,ni) &
2754  + mode/ray*jsolh_anal(6)*h_mesh%gauss%ww(ni,l) &
2755  + rhs_hl(1)*h_mesh%gauss%ww(ni,l))
2756 
2757  src_h(i,2) = src_h(i,2)+ ray_rjl &
2758  *(jsolh_anal(4)*dwh(2,ni) &
2759  - mode/ray*jsolh_anal(5)*h_mesh%gauss%ww(ni,l) &
2760  + rhs_hl(2)*h_mesh%gauss%ww(ni,l))
2761 
2762  !--------Composante theta------
2763  src_h(i,3) = src_h(i,3)+ ray_rjl &
2764  * (-jsolh_anal(1)*dwh(2,ni) &
2765  + 1/ray*jsolh_anal(5)*(ray*dwh(1,ni) + h_mesh%gauss%ww(ni,l)) &
2766  + rhs_hl(3)*h_mesh%gauss%ww(ni,l))
2767 
2768  src_h(i,4) = src_h(i,4)+ ray_rjl &
2769  * (-jsolh_anal(2)*dwh(2,ni) &
2770  + 1/ray*jsolh_anal(6)*(ray*dwh(1,ni) + h_mesh%gauss%ww(ni,l)) &
2771  + rhs_hl(4)*h_mesh%gauss%ww(ni,l))
2772 
2773  !--------Composante z------
2774  src_h(i,5) = src_h(i,5)+ ray_rjl* &
2775  (-mode/ray*jsolh_anal(2)*h_mesh%gauss%ww(ni,l) &
2776  - jsolh_anal(3)*dwh(1,ni) &
2777  + rhs_hl(5)*h_mesh%gauss%ww(ni,l))
2778 
2779  src_h(i,6) = src_h(i,6)+ ray_rjl* &
2780  (mode/ray*jsolh_anal(1)*h_mesh%gauss%ww(ni,l) &
2781  - jsolh_anal(4)*dwh(1,ni) &
2782  + rhs_hl(6)*h_mesh%gauss%ww(ni,l))
2783  ENDDO
2784 
2785  END DO
2786  END DO
2787  !tps = user_time(dummy)- tps
2788  !WRITE(*,*) ' Temps in courant boucle me H', tps
2789  !tps = user_time(dummy)
2790 
2791  ! We integrate by parts this term
2792  ! JLG + FL, FEB 10, 2010
2793  !DO m = 1, phi_mesh%me
2794  ! src_phil=0
2795  ! DO l = 1, phi_mesh%gauss%l_G
2796  ! dwphi = phi_mesh%gauss%dw(:,:,l,m)
2797  ! !===Compute radius of Gauss point
2798  ! rhs_dphil=0
2799  ! rhs_phil=0
2800  ! ray = 0
2801  ! DO ni = 1, phi_mesh%gauss%n_w; i = phi_mesh%jj(ni,m)
2802  ! ray = ray + phi_mesh%rr(1,i)*phi_mesh%gauss%ww(ni,l)
2803  ! rhs_phil(:) = rhs_phil(:) + rhs_phi(i,:)*phi_mesh%gauss%ww(ni,l)
2804  ! DO k =1 ,2
2805  ! rhs_dphil(:,k) = rhs_dphil(:,k) + rhs_phi(i,:)*dwphi(k,ni)
2806  ! END DO
2807  ! END DO
2808  ! ray_rjl = phi_mesh%gauss%rj(l,m)*ray
2809  ! moderay2 = (mode/ray)**2
2810 
2811  ! DO ni = 1, phi_mesh%gauss%n_w
2812 
2813  ! src_phil(1,ni) = src_phil(1,ni) + ray_rjl* &
2814  ! (rhs_dphil(1,1)*dwphi(1,ni) + &
2815  ! moderay2*rhs_phil(1)*phi_mesh%gauss%ww(ni,l) + &
2816  ! rhs_dphil(1,2)*dwphi(2,ni))
2817 
2818  ! src_phil(2,ni) = src_phil(2,ni) + ray_rjl* &
2819  ! (rhs_dphil(2,1)*dwphi(1,ni) + &
2820  ! moderay2*rhs_phil(2)*phi_mesh%gauss%ww(ni,l) + &
2821  ! rhs_dphil(2,2)*dwphi(2,ni))
2822  ! END DO
2823 
2824  ! END DO
2825  ! DO ni = 1, phi_mesh%gauss%n_w
2826  ! i = phi_mesh%jj(ni,m)
2827  ! src_phi(i,:) = src_phi(i,:) + src_phil(:,ni)
2828  ! END DO
2829  !END DO
2830  ! End integration by parts
2831  ! JLG + FL, FEB 10, 2010
2832 
2833 
2834  !==interface_H_phi
2835  !forcage sur l'interface_H_phi
2836  !attention on comprime le calcul sur les points de Gauss et integration !!
2837  !j/sigma*(b x nc + grad(phi) x nv)
2838 
2839  CALL gauss(phi_mesh)
2840 
2841  n_ws1 = h_mesh%gauss%n_ws
2842  n_ws2 = phi_mesh%gauss%n_ws
2843  n_w2 = phi_mesh%gauss%n_w
2844 
2845  h_bloc_size = h_mesh%np
2846 
2847  IF (interface_h_phi%mes /=0) THEN ! Ajout du test pour les grands nb de domaines
2848  IF (h_mesh%gauss%n_ws == n_ws) THEN
2849  w_cs = wws
2850  ELSE
2851  DO ls = 1, l_gs
2852  w_cs(1,ls)= wws(1,ls)+0.5*wws(3,ls)
2853  w_cs(2,ls)= wws(2,ls)+0.5*wws(3,ls)
2854  w_cs(3,ls)=0
2855  ENDDO
2856  END IF
2857  END IF
2858 
2859 
2860  !WRITE(*,*) ' Courant: init gauss'
2861  DO ms = 1, interface_h_phi%mes
2862 
2863  ms2 = interface_h_phi%mesh2(ms)
2864  ms1 = interface_h_phi%mesh1(ms)
2865  m = phi_mesh%neighs(ms2)
2866  m1 = h_mesh%neighs(ms1)
2867  mesh_id1 = h_mesh%i_d(m1)
2868  DO ls = 1,l_gs
2869  !Feb 9 2007, muhl
2870  muhl=sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
2871  !Feb 9 2007, muhl
2872  DO k=1, 6
2873  h_ext_l(k) = sum(h_ext(interface_h_phi%jjs1(1:n_ws1,ms),k)*w_cs(1:n_ws1,ls))
2874  j_over_sigma_gauss(k) = sum(j_over_sigma(interface_h_phi%jjs1(1:n_ws1,ms),k)*w_cs(1:n_ws1,ls))
2875  END DO
2876 
2877  !===Compute radius of Gauss point
2878  ray = 0
2879  DO ni = 1, n_ws2; i = phi_mesh%jjs(ni,interface_h_phi%mesh2(ms))
2880  ray = ray + phi_mesh%rr(1,i)* wws(ni,ls)
2881  END DO
2882 
2883  gaussp = 0.d0
2884  DO ni=1, n_ws2
2885  i=phi_mesh%jjs(ni,ms2)
2886  gaussp = gaussp + phi_mesh%rr(:,i)*phi_mesh%gauss%wws(ni,ls)
2887  ENDDO
2888 
2889  DO k=1, 6
2890  jsolh_anal(k) = jexact_gauss(k, gaussp, mode, mu_phi ,sigma(m1), &
2891  muhl, time, mesh_id1, h_ext_l)/sigma(m1) &
2892  + muhl * sum(nl(h_mesh%jjs(1:n_ws1,ms1),k)*w_cs(1:n_ws1,ls))
2893  ENDDO
2894 !!$! TEST DEBUG : to do before using H with phi
2895 !!$ IF (inputs%if_level_set.AND.inputs%variation_sigma_fluid) THEN
2896 !!$ DO k = 1, 6
2897 !!$ JsolH_anal(k) = J_over_sigma_gauss(k) + sigma_curl(index,k) &
2898 !!$ + SUM(NL(H_mesh%jjs(1:n_ws1,ms1),k)*w_cs(1:n_ws1,ls))
2899 !!$ END DO
2900 !!$ ELSE
2901 !!$ DO k = 1, 6
2902 !!$ JsolH_anal(k) = Jexact_gauss(k, gaussp, mode, mu_phi ,sigma(m1), &
2903 !!$ muhl, time, mesh_id1, H_ext_l)/sigma(m1) &
2904 !!$ + SUM(NL(H_mesh%jjs(1:n_ws1,ms1),k)*w_cs(1:n_ws1,ls))
2905 !!$ END DO
2906 !!$ END IF
2907 !!$! TEST DEBUG
2908 
2909  !---------forcage pour H
2910 
2911  DO ni=1, n_ws1
2912  i = interface_h_phi%jjs1(ni,ms)
2913  src_h(i,1) = src_h(i,1)+rjs(ls,ms2)*ray*( &
2914  -jsolh_anal(3)*w_cs(ni,ls)*(-rnorms(2,ls,ms2)))
2915 
2916  src_h(i,2) = src_h(i,2)+rjs(ls,ms2)*ray*( &
2917  -jsolh_anal(4)*w_cs(ni,ls)*(-rnorms(2,ls,ms2)))
2918 
2919  src_h(i,3) = src_h(i,3)+rjs(ls,ms2)*ray*( &
2920  jsolh_anal(1)*w_cs(ni,ls)*(-rnorms(2,ls,ms2)) &
2921  -jsolh_anal(5)*w_cs(ni,ls)*(-rnorms(1,ls,ms2)))
2922 
2923  src_h(i,4) = src_h(i,4)+rjs(ls,ms2)*ray*( &
2924  jsolh_anal(2)*w_cs(ni,ls)*(-rnorms(2,ls,ms2)) &
2925  -jsolh_anal(6)*w_cs(ni,ls)*(-rnorms(1,ls,ms2)))
2926 
2927  src_h(i,5) = src_h(i,5)+rjs(ls,ms2)*ray*( &
2928  jsolh_anal(3)*w_cs(ni,ls)*(-rnorms(1,ls,ms2)))
2929 
2930  src_h(i,6) = src_h(i,6)+rjs(ls,ms2)*ray*( &
2931  jsolh_anal(4)*w_cs(ni,ls)*(-rnorms(1,ls,ms2)))
2932  ENDDO
2933 
2934  !---------forcage pour phi
2935  !terme sans derivee de phi
2936  DO ni=1,n_ws2
2937  i = interface_h_phi%jjs2(ni,ms)
2938  !attention si on force sur l'axe, il faut retirer les 1/ray
2939  !There was a BUG here. There was w_cs instead of wws
2940  src_phi(i,1) = src_phi(i,1)+rjs(ls,ms2)*( &
2941  - mode*jsolh_anal(2)*wws(ni,ls) * rnorms(2,ls,ms2) &
2942  + mode*jsolh_anal(6)*wws(ni,ls) * rnorms(1,ls,ms2))
2943 
2944  src_phi(i,2) = src_phi(i,2)+rjs(ls,ms2)*( &
2945  + mode*jsolh_anal(1)*wws(ni,ls) * rnorms(2,ls,ms2) &
2946  - mode*jsolh_anal(5)*wws(ni,ls) * rnorms(1,ls,ms2))
2947 
2948  ENDDO
2949 
2950  !terme avec derivee de phi
2951  DO ni=1,n_w2
2952  i = phi_mesh%jj(ni,m)
2953  src_phi(i,1) = src_phi(i,1)+rjs(ls,ms2)*ray*( &
2954  + jsolh_anal(3) *(dw_s(2,ni,ls,ms2) * rnorms(1,ls,ms2)&
2955  -dw_s(1,ni,ls,ms2) * rnorms(2,ls,ms2)))
2956 
2957  src_phi(i,2) = src_phi(i,2)+rjs(ls,ms2)*ray*( &
2958  + jsolh_anal(4)*(dw_s(2,ni,ls,ms2) * rnorms(1,ls,ms2)&
2959  -dw_s(1,ni,ls,ms2) * rnorms(2,ls,ms2)))
2960 
2961  ENDDO
2962 
2963  ! Integration by parts of int(GRAD rhs_phi Grad psi)
2964  rhs_hl = 0.d0
2965  DO ni=1, n_ws1
2966  i = interface_h_phi%jjs1(ni,ms)
2967  rhs_hl(:) = rhs_hl(:) + rhs_h(i,:)*w_cs(ni,ls)
2968  ENDDO
2969 
2970  DO ni=1, n_ws2
2971  i = interface_h_phi%jjs2(ni,ms)
2972  src_phi(i,1) = src_phi(i,1)+rjs(ls,ms2)*ray*wws(ni,ls)*( &
2973  rhs_hl(1)*rnorms(1,ls,ms2) + rhs_hl(5)*rnorms(2,ls,ms2))
2974  src_phi(i,2) = src_phi(i,2)+rjs(ls,ms2)*ray*wws(ni,ls)*( &
2975  rhs_hl(2)*rnorms(1,ls,ms2) + rhs_hl(6)*rnorms(2,ls,ms2))
2976  END DO
2977  ! End integration by parts of int(GRAD rhs_phi Grad psi)
2978  END DO
2979  END DO
2980  !tps = user_time(dummy)- tps
2981  !WRITE(*,*) ' Courant: init interface_H_phi'
2982 
2983  IF (h_mesh%np /= 0) THEN
2984 !!$ ALLOCATE(idxn(H_mesh%np))
2985  idxn_h = la_h%loc_to_glob(1,:)-1
2986  CALL vecsetvalues(vb_1, h_mesh%np, idxn_h, src_h(:,1), add_values, ierr)
2987  CALL vecsetvalues(vb_2, h_mesh%np, idxn_h, src_h(:,2), add_values, ierr)
2988  idxn_h = la_h%loc_to_glob(2,:)-1
2989  CALL vecsetvalues(vb_1, h_mesh%np, idxn_h, src_h(:,4), add_values, ierr)
2990  CALL vecsetvalues(vb_2, h_mesh%np, idxn_h, src_h(:,3), add_values, ierr)
2991  idxn_h = la_h%loc_to_glob(3,:)-1
2992  CALL vecsetvalues(vb_1, h_mesh%np, idxn_h, src_h(:,5), add_values, ierr)
2993  CALL vecsetvalues(vb_2, h_mesh%np, idxn_h, src_h(:,6), add_values, ierr)
2994 !!$ DEALLOCATE(idxn)
2995  END IF
2996  IF (phi_mesh%np /=0) THEN
2997 !!$ ALLOCATE(idxn(phi_mesh%np))
2998  idxn_phi = la_phi%loc_to_glob(1,:)-1
2999  CALL vecsetvalues(vb_1, phi_mesh%np, idxn_phi, src_phi(:,1), add_values, ierr)
3000  CALL vecsetvalues(vb_2, phi_mesh%np, idxn_phi, src_phi(:,2), add_values, ierr)
3001 !!$ DEALLOCATE(idxn)
3002  END IF
3003 
3004  CALL vecassemblybegin(vb_1,ierr)
3005  CALL vecassemblyend(vb_1,ierr)
3006  CALL vecassemblybegin(vb_2,ierr)
3007  CALL vecassemblyend(vb_2,ierr)
3008 
3009 !!$ IF (H_mesh%me /=0) THEN
3010 !!$ DEALLOCATE(nl)
3011 !!$ END IF
3012  !DEALLOCATE(src_H, src_phi)
3013 
3014  END SUBROUTINE courant_int_by_parts
3015 
3016  SUBROUTINE surf_int(H_mesh,phi_mesh, interface_H_phi, interface_H_mu, &
3017  list_dirichlet_sides_h, sigma,mu_phi, mu_h_field,time,mode, &
3018  la_h, la_phi, vb_1, vb_2, r_fourier, index_fourier)
3019  !calcul du forcage a la frontiere exterieure
3020  USE my_util
3021  USE def_type_mesh
3022  USE boundary
3023  USE input_data
3024  IMPLICIT NONE
3025  TYPE(mesh_type), INTENT(IN) :: h_mesh, phi_mesh
3026  TYPE(interface_type), INTENT(IN) :: interface_h_phi, interface_h_mu
3027  INTEGER, DIMENSION(:), INTENT(IN) :: list_dirichlet_sides_h
3028  REAL(KIND=8), INTENT(IN) :: mu_phi, time
3029  REAL(KIND=8),DIMENSION(H_mesh%me),INTENT(IN):: sigma
3030  REAL(KIND=8),DIMENSION(:), INTENT(IN) :: mu_h_field
3031  INTEGER, INTENT(IN) :: mode
3032  REAL(KIND=8), OPTIONAL :: r_fourier
3033  INTEGER, OPTIONAL :: index_fourier
3034  REAL(KIND=8), DIMENSION(H_mesh%np, 6) :: src_h
3035  REAL(KIND=8), DIMENSION(phi_mesh%np, 2) :: src_phi
3036 ! CN possible faute POINTER ?
3037  !REAL(KIND=8), DIMENSION(:,:), POINTER :: src_H, src_phi
3038  !REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: src_H, src_phi
3039  REAL(KIND=8) :: ray, muhl
3040  INTEGER :: ms, ls, ns, i, k, m, n, ni
3041  REAL(KIND=8), DIMENSION(2) :: gaussp
3042  REAL(KIND=8), DIMENSION(6) :: esolh_anal, esolphi_anal
3043 !!$ FL+CN 22/03/2013
3044 !!$ INTEGER, DIMENSION(:), ALLOCATABLE :: idxn
3045  INTEGER, DIMENSION(H_mesh%np) :: idxn_h
3046  INTEGER, DIMENSION(phi_mesh%np) :: idxn_phi
3047 !!$ FL+CN 22/03/2013
3048  TYPE(petsc_csr_la) :: la_h, la_phi
3049 
3050  INTEGER :: count
3051  INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: neumann_bdy_h_sides
3052  INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: neumann_bdy_phi_sides
3053  LOGICAL, SAVE :: once_neumann_bdy=.true.
3054  LOGICAL, DIMENSION(:), ALLOCATABLE :: virgin1
3055  LOGICAL, DIMENSION(:), ALLOCATABLE :: virgin2
3056 
3057 #include "petsc/finclude/petsc.h"
3058  petscerrorcode :: ierr
3059  vec :: vb_1, vb_2
3060 
3061  !ALLOCATE(src_H(H_mesh%np,6), src_phi(phi_mesh%np,2))
3062 
3063  !===Compute sides that are on Neumann boundary:
3064  !===Neumann_bdy_H_sides, Neumann_bdy_phi_sides
3065  IF (once_neumann_bdy) THEN
3066  once_neumann_bdy = .false.
3067  ALLOCATE(virgin1(h_mesh%mes))
3068  ALLOCATE(virgin2(phi_mesh%mes))
3069  virgin1=.true.
3070  virgin2=.true.
3071  IF (interface_h_phi%mes/=0) THEN
3072  virgin1(interface_h_phi%mesh1) = .false.
3073  virgin2(interface_h_phi%mesh2) = .false.
3074  END IF
3075  IF (interface_h_mu%mes/=0) THEN
3076  virgin1(interface_h_mu%mesh1) = .false.
3077  virgin1(interface_h_mu%mesh2) = .false.
3078  END IF
3079  !===Create Neumann_bdy_H_sides
3080  count = 0
3081  DO ms = 1, h_mesh%mes
3082  IF (maxval(abs(h_mesh%rr(1,h_mesh%jjs(:,ms)))).LT.1d-12) cycle ! No Neumann BC on the z-axis
3083  IF (.NOT.virgin1(ms)) cycle ! No Neumann BC on H_mu interface
3084  IF(minval(abs(h_mesh%sides(ms)-list_dirichlet_sides_h))==0) cycle ! Dirichlet boundary
3085  count = count + 1
3086  END DO
3087  ALLOCATE(neumann_bdy_h_sides(count))
3088  count = 0
3089  DO ms = 1, h_mesh%mes
3090  IF (maxval(abs(h_mesh%rr(1,h_mesh%jjs(:,ms)))).LT.1d-12) cycle
3091  IF (.NOT.virgin1(ms)) cycle
3092  IF(minval(abs(h_mesh%sides(ms)-list_dirichlet_sides_h))==0) cycle
3093  count = count + 1
3094  neumann_bdy_h_sides(count) = ms
3095  END DO
3096  !===Create Neumann_bdy_phi_sides
3097  count = 0
3098  DO ms = 1, phi_mesh%mes
3099  IF (present(index_fourier)) THEN
3100  IF (phi_mesh%sides(ms)==index_fourier) cycle ! No Neumann BC on Fourier boundary
3101  END IF
3102  IF (.NOT.virgin2(ms)) cycle ! No Neumann BC on H_phi interface
3103  IF (maxval(abs(phi_mesh%rr(1,phi_mesh%jjs(:,ms)))).LT.1d-12) cycle ! No Neumann BC on the z-axis
3104  IF (minval(abs(phi_mesh%sides(ms)-inputs%phi_list_dirichlet_sides))==0) cycle ! Dirichlet boundary
3105  count = count + 1
3106  END DO
3107  ALLOCATE(neumann_bdy_phi_sides(count))
3108  count = 0
3109  DO ms = 1, phi_mesh%mes
3110  IF (present(index_fourier)) THEN
3111  IF (phi_mesh%sides(ms)==index_fourier) cycle
3112  END IF
3113  IF (.NOT.virgin2(ms)) cycle
3114  IF (maxval(abs(phi_mesh%rr(1,phi_mesh%jjs(:,ms)))).LT.1d-12) cycle
3115  IF (minval(abs(phi_mesh%sides(ms)-inputs%phi_list_dirichlet_sides))==0) cycle ! Dirichlet boundary
3116  count = count + 1
3117  neumann_bdy_phi_sides(count) = ms
3118  END DO
3119  DEALLOCATE(virgin1)
3120  DEALLOCATE(virgin2)
3121  END IF
3122 
3123  src_h = 0.d0
3124  src_phi = 0.d0
3125 
3126  DO count = 1, SIZE(neumann_bdy_h_sides)
3127  ms = neumann_bdy_h_sides(count)
3128  m = h_mesh%neighs(ms)
3129  DO ls = 1, h_mesh%gauss%l_Gs
3130  !Feb 8 2007, mmuhl
3131  muhl = sum(mu_h_field(h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
3132  !Feb 8 2007, mmuhl
3133 
3134  !===Compute radius of Gauss point
3135  ray = 0
3136  DO ni = 1, h_mesh%gauss%n_ws; i = h_mesh%jjs(ni,ms)
3137  ray = ray + h_mesh%rr(1,i)* h_mesh%gauss%wws(ni,ls)
3138  END DO
3139 
3140  IF (ray.LT.1.d-15) cycle !ATTENTION Axe
3141 
3142  gaussp = 0.d0
3143  DO ns=1, h_mesh%gauss%n_ws
3144  i=h_mesh%jjs(ns,ms)
3145  gaussp = gaussp + h_mesh%rr(:,i)*h_mesh%gauss%wws(ns,ls)
3146  ENDDO
3147 
3148  DO k=1, 6
3149  esolh_anal(k) = eexact_gauss(k,gaussp,mode,mu_phi,sigma(m),muhl, time)
3150  ENDDO
3151 
3152  !forcage pour la frontiere de H
3153  ! - E.(b x nc)
3154 
3155  DO ns=1, h_mesh%gauss%n_ws
3156  i = h_mesh%jjs(ns,ms)
3157  src_h(i,1) = src_h(i,1)-h_mesh%gauss%rjs(ls,ms)*ray*( &
3158  -esolh_anal(3)*h_mesh%gauss%wws(ns,ls)* &
3159  (h_mesh%gauss%rnorms(2,ls,ms)))
3160 
3161  src_h(i,2) = src_h(i,2)-h_mesh%gauss%rjs(ls,ms)*ray*( &
3162  -esolh_anal(4)*h_mesh%gauss%wws(ns,ls)* &
3163  (h_mesh%gauss%rnorms(2,ls,ms)))
3164 
3165  src_h(i,3) = src_h(i,3)-h_mesh%gauss%rjs(ls,ms)*ray*( &
3166  esolh_anal(1)*h_mesh%gauss%wws(ns,ls)* &
3167  (h_mesh%gauss%rnorms(2,ls,ms)) - &
3168  esolh_anal(5)*h_mesh%gauss%wws(ns,ls) * &
3169  (h_mesh%gauss%rnorms(1,ls,ms)))
3170 
3171  src_h(i,4) = src_h(i,4)-h_mesh%gauss%rjs(ls,ms)*ray*( &
3172  esolh_anal(2)*h_mesh%gauss%wws(ns,ls) * &
3173  (h_mesh%gauss%rnorms(2,ls,ms)) - &
3174  esolh_anal(6)*h_mesh%gauss%wws(ns,ls) * &
3175  (h_mesh%gauss%rnorms(1,ls,ms)))
3176 
3177  src_h(i,5) = src_h(i,5)-h_mesh%gauss%rjs(ls,ms)*ray*( &
3178  esolh_anal(3)*h_mesh%gauss%wws(ns,ls)* &
3179  (h_mesh%gauss%rnorms(1,ls,ms)))
3180 
3181  src_h(i,6) = src_h(i,6)-h_mesh%gauss%rjs(ls,ms)*ray*( &
3182  esolh_anal(4)*h_mesh%gauss%wws(ns,ls) * &
3183  (h_mesh%gauss%rnorms(1,ls,ms)))
3184 
3185  ENDDO
3186 
3187  ENDDO
3188  ENDDO
3189 
3190  !===Neumann boundary phi_mesh
3191  DO count = 1, SIZE(neumann_bdy_phi_sides)
3192  ms = neumann_bdy_phi_sides(count)
3193  m = phi_mesh%neighs(ms)
3194  DO ls = 1, phi_mesh%gauss%l_Gs
3195  !===Compute radius of Gauss point
3196  ray = 0
3197  DO ni = 1, phi_mesh%gauss%n_ws; i = phi_mesh%jjs(ni,ms)
3198  ray = ray + phi_mesh%rr(1,i)* phi_mesh%gauss%wws(ni,ls)
3199  END DO
3200 
3201  gaussp = 0.d0
3202  DO ns=1, phi_mesh%gauss%n_ws
3203  i=phi_mesh%jjs(ns,ms)
3204  gaussp = gaussp + phi_mesh%rr(:,i)*phi_mesh%gauss%wws(ns,ls)
3205  ENDDO
3206 
3207  DO k=1, 6
3208  !===Here sigma and mu_H_field should not intervene
3209  !===I put boggus values for sigma and mu_H_field, Feb 8 2007, Jean-Luc Guermond
3210  esolphi_anal(k) = eexact_gauss(k,gaussp,mode,mu_phi,sigma(1),mu_h_field(1),time)
3211  ENDDO
3212 
3213  !===Nemnann forcing for phi in rhs: - E.(grad(phi) x nv)
3214  DO ns=1, phi_mesh%gauss%n_ws
3215  i = phi_mesh%jjs(ns,ms)
3216  DO n = 1, phi_mesh%gauss%n_w
3217  IF (phi_mesh%jj(n,m) == i) EXIT
3218  END DO
3219  !===There should not be any Neumann forcing on z-axis (1/ray would be infinite)
3220  src_phi(i,1) = src_phi(i,1)-phi_mesh%gauss%rjs(ls,ms)*ray*( &
3221  +esolphi_anal(3)*(phi_mesh%gauss%dw_s(2,n,ls,ms)*phi_mesh%gauss%rnorms(1,ls,ms) &
3222  -phi_mesh%gauss%dw_s(1,n,ls,ms)*phi_mesh%gauss%rnorms(2,ls,ms))) &
3223  -phi_mesh%gauss%rjs(ls,ms)*(&
3224  -mode*esolphi_anal(2)*phi_mesh%gauss%wws(ns,ls)*phi_mesh%gauss%rnorms(2,ls,ms) &
3225  +mode*esolphi_anal(6)*phi_mesh%gauss%wws(ns,ls)*phi_mesh%gauss%rnorms(1,ls,ms))
3226 
3227  src_phi(i,2) = src_phi(i,2)-phi_mesh%gauss%rjs(ls,ms)*ray*( &
3228  +esolphi_anal(4)*(phi_mesh%gauss%dw_s(2,n,ls,ms)*phi_mesh%gauss%rnorms(1,ls,ms) &
3229  -phi_mesh%gauss%dw_s(1,n,ls,ms)*phi_mesh%gauss%rnorms(2,ls,ms))) &
3230  -phi_mesh%gauss%rjs(ls,ms)*(&
3231  mode*esolphi_anal(1)*phi_mesh%gauss%wws(ns,ls)*phi_mesh%gauss%rnorms(2,ls,ms) &
3232  -mode*esolphi_anal(5)*phi_mesh%gauss%wws(ns,ls)*phi_mesh%gauss%rnorms(1,ls,ms))
3233  ENDDO
3234  ENDDO
3235  ENDDO
3236  !===ATTENTION ATTENTION ATTENTION ATTENTION ATTENTION ATTENTION ATTENTION ATTENTION
3237  !JLG, FL, FEB, 10, 2010
3238  !We assume that integral int_Gammav n.GRAD (4phi_n-phi_(n-1))/(2dt) psi dsigma = 0
3239  !JLG, FL, FEB, 10, 2010
3240  !JLG, FL, May, 28, 2009
3241  !We assume that integral int_Gammav (Hinfty . n) psi ds = 0!
3242  !JLG, FL, May, 28, 2009
3243  !===ATTENTION ATTENTION ATTENTION ATTENTION ATTENTION ATTENTION ATTENTION ATTENTION
3244 
3245  !=========================================================
3246  !--- Artificial boundary condition: d(phi_t)/dR + (1/R)*phi_t = 0
3247  !=========================================================
3248 
3249  IF (present(r_fourier)) THEN
3250  IF (r_fourier.GE.0.d0) CALL error_petsc('maxwell_update_time_with_H: R_fourier should be -1')
3251  END IF
3252  !IF (.NOT.present(index_fourier) .OR. .NOT.present(R_fourier)) RETURN
3253  !IF (R_fourier.le.0.d0) RETURN
3254  !DO ms = 1, phi_mesh%mes
3255  ! IF (phi_mesh%sides(ms) /= index_fourier) CYCLE ! Not on the artificial boundary
3256 
3257  ! DO ls = 1, phi_mesh%gauss%l_Gs
3258 
3259  !===Compute radius of Gauss point
3260  ! ray = SUM(phi_mesh%rr(1,phi_mesh%jjs(:,ms))* phi_mesh%gauss%wws(:,ls))
3261  ! x = phi_mesh%gauss%rjs(ls,ms)*ray/R_fourier
3262  ! y1 = x* SUM(rhs(phi_mesh%jjs(:,ms),1)* phi_mesh%gauss%wws(:,ls))
3263  ! y2 = x* SUM(rhs(phi_mesh%jjs(:,ms),2)* phi_mesh%gauss%wws(:,ls))
3264  ! DO ns =1, phi_mesh%gauss%n_ws
3265  ! src_phi(1,phi_mesh%jjs(ns,ms)) = src_phi(1,phi_mesh%jjs(ns,ms)) + &
3266  ! y1*phi_mesh%gauss%wws(ns,ls)
3267  ! src_phi(2,phi_mesh%jjs(ns,ms)) = src_phi(2,phi_mesh%jjs(ns,ms)) + &
3268  ! y2*phi_mesh%gauss%wws(ns,ls)
3269  ! ENDDO
3270  !
3271  ! ENDDO
3272  !END DO
3273  IF (h_mesh%mes/=0) THEN
3274 !!$ ALLOCATE(idxn(H_mesh%np))
3275  idxn_h = la_h%loc_to_glob(1,:)-1
3276  CALL vecsetvalues(vb_1, h_mesh%np, idxn_h, src_h(:,1), add_values, ierr)
3277  CALL vecsetvalues(vb_2, h_mesh%np, idxn_h, src_h(:,2), add_values, ierr)
3278  idxn_h = la_h%loc_to_glob(2,:)-1
3279  CALL vecsetvalues(vb_1, h_mesh%np, idxn_h, src_h(:,4), add_values, ierr)
3280  CALL vecsetvalues(vb_2, h_mesh%np, idxn_h, src_h(:,3), add_values, ierr)
3281  idxn_h = la_h%loc_to_glob(3,:)-1
3282  CALL vecsetvalues(vb_1, h_mesh%np, idxn_h, src_h(:,5), add_values, ierr)
3283  CALL vecsetvalues(vb_2, h_mesh%np, idxn_h, src_h(:,6), add_values, ierr)
3284 !!$ DEALLOCATE(idxn)
3285  END IF
3286  IF (phi_mesh%mes/=0) THEN
3287 !!$ ALLOCATE(idxn(phi_mesh%np))
3288  idxn_phi = la_phi%loc_to_glob(1,:)-1
3289  CALL vecsetvalues(vb_1, phi_mesh%np, idxn_phi, src_phi(:,1), add_values, ierr)
3290  CALL vecsetvalues(vb_2, phi_mesh%np, idxn_phi, src_phi(:,2), add_values, ierr)
3291 !!$ DEALLOCATE(idxn)
3292  END IF
3293 
3294  CALL vecassemblybegin(vb_1,ierr)
3295  CALL vecassemblyend(vb_1,ierr)
3296  CALL vecassemblybegin(vb_2,ierr)
3297  CALL vecassemblyend(vb_2,ierr)
3298 
3299  !DEALLOCATE(src_H,src_phi)
3300 
3301  END SUBROUTINE surf_int
3302 
3303  SUBROUTINE mat_maxwell_mu(H_mesh, interface_H_mu, mode, stab, &
3304  mu_h_field, sigma, la_h, h_p_phi_mat1, h_p_phi_mat2)
3305  USE def_type_mesh
3306  USE gauss_points
3307  USE my_util
3308  IMPLICIT NONE
3309  TYPE(mesh_type), INTENT(IN) :: h_mesh
3310  TYPE(interface_type), INTENT(IN) :: interface_h_mu
3311  INTEGER, INTENT(IN) :: mode
3312  REAL(KIND=8), DIMENSION(3), INTENT(IN) :: stab
3313  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: sigma, mu_h_field
3314  INTEGER :: ms, ls, ni, nj, i, j, &
3315  n_ws1, n_ws2, n_w2, n_w1, m1, m2, ki, kj,ib,jb, ms1, ms2
3316  REAL(KIND=8) :: x, y, z, norm, hm1
3317  REAL(KIND=8) :: ray, stab_colle_h_mu
3318  LOGICAL :: mark=.false.
3319  REAL(KIND=8), DIMENSION(9,SIZE(H_mesh%jj,1),SIZE(H_mesh%jj,1),2,2) :: hsij, gsij
3320 
3321  ! MATRICES POUR LES TERMES DE BORDS Hsij et Gsij
3322  !=================================================
3323  ! (--------------------------------------------------)
3324  ! ( Hsij(1) +G | GSij(2) | Hsij(4) +G )
3325  ! ( Hsij(1) +G | GSij(3) | Hsij(4) +G )
3326  ! (--------------------------------------------------)
3327  ! ( Hsij(2) | Hsij(5) +G | Hsij(8) )
3328  ! ( Hsij(3) | Hsij(5) +G | Hsij(9) )
3329  ! (--------------------------------------------------)
3330  ! ( Hsij(7) +G | GSij(8) | Hsij(6) +G )
3331  ! ( Hsij(7) +G | GSij(9) | Hsij(6) +G )
3332  ! (==================================================)
3333 !!$ FL+CN 22/03/2013
3334  REAL(KIND=8), DIMENSION(H_mesh%gauss%n_ws,H_mesh%gauss%l_Gs) :: w_cs
3335  REAL(KIND=8), DIMENSION(2, H_mesh%gauss%n_w, H_mesh%gauss%l_Gs, H_mesh%mes) :: dw_cs
3336  REAL(KIND=8), DIMENSION(2, H_mesh%gauss%n_w) :: dwsi,dwsj
3337  REAL(KIND=8), DIMENSION(2,H_mesh%gauss%l_Gs) :: gauss1, gauss2
3338 !!$ REAL(KIND=8), DIMENSION(:,:) , ALLOCATABLE :: w_cs
3339 !!$ REAL(KIND=8), DIMENSION(:,:,:,:), ALLOCATABLE :: dw_cs
3340 !!$ REAL(KIND=8), DIMENSION(:,:) , ALLOCATABLE :: dwsi,dwsj
3341 !!$ REAL(KIND=8), DIMENSION(:,:) , ALLOCATABLE :: gauss1, gauss2
3342 !!$ FL+CN 22/03/2013
3343  REAL(KIND=8), DIMENSION(2) :: normi, normj
3344  REAL(KIND=8), DIMENSION(SIZE(H_mesh%jjs,1)) :: wwsi, wwsj
3345  INTEGER :: n_wsi, n_wsj, ci, cj, n_wi, n_wj
3346 
3347  INTEGER :: ls1, ls2
3348  REAL(KIND=8) :: ref, diff, mu_h, muhl1, muhl2, muhi, muhj, sigmai, sigmaj
3349  ! June 14 2008
3350  REAL(KIND=8) :: c_sym =.0d0 ! (c_sym=1.d0 symmetrizes the bilinear form)
3351  REAL(KIND=8) :: wwiwwj, normt, stab_div
3352  ! June 14 2008
3353 !!$ FL +CN 22/03/2013
3354 !!$ REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: mat_loc1, mat_loc2
3355 !!$ INTEGER , DIMENSION(:), ALLOCATABLE :: idxn, jdxn
3356  REAL(KIND=8), DIMENSION(6*H_mesh%gauss%n_w,6*H_mesh%gauss%n_w) :: mat_loc1, mat_loc2
3357  INTEGER , DIMENSION(6*H_mesh%gauss%n_w) :: idxn, jdxn
3358 !!$ FL +CN 22/03/2013
3359  TYPE(petsc_csr_la) :: la_h
3360  INTEGER :: ix, jx
3361 
3362 #include "petsc/finclude/petsc.h"
3363  mat :: h_p_phi_mat1, h_p_phi_mat2
3364  petscerrorcode :: ierr
3365 
3366  ! June 2009, JLG, CN, normalization
3367  stab_colle_h_mu = stab(3)
3368  stab_div = stab(1)
3369  ! June 2009, JLG, CN
3370 
3371  !**********************************************************************************
3372  !--------------------TERMS ON SIGMA_MU-------------------------------
3373  !**********************************************************************************
3374 
3375  !WRITE(*,*) 'Assembling interface_H_mu'
3376  CALL gauss(h_mesh)
3377  n_ws1 = h_mesh%gauss%n_ws
3378  n_ws2 = h_mesh%gauss%n_ws
3379  n_w1 = h_mesh%gauss%n_w
3380  n_w2 = h_mesh%gauss%n_w
3381 
3382 !!$ ALLOCATE(w_cs(n_ws1,l_Gs))
3383 !!$ ALLOCATE(dw_cs(2, n_w1, l_Gs, H_mesh%mes))
3384 !!$ ALLOCATE(dwsi(2, n_w1),dwsj(2, n_w2))
3385 !!$ ALLOCATE(gauss1(2,l_Gs),gauss2(2,l_Gs))
3386 
3387 !!$ ALLOCATE(mat_loc1(6*n_w1,6*n_w2))
3388 !!$ ALLOCATE(mat_loc2(6*n_w1,6*n_w2))
3389 !!$ ALLOCATE(idxn(6*n_w1))
3390 !!$ ALLOCATE(jdxn(6*n_w2))
3391 
3392  DO ms = 1, interface_h_mu%mes
3393 
3394  ms2 = interface_h_mu%mesh2(ms)
3395  m2 = h_mesh%neighs(ms2)
3396  ms1 = interface_h_mu%mesh1(ms)
3397  m1 = h_mesh%neighs(ms1)
3398 
3399  ref = 1.d-8+sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - h_mesh%rr(:,h_mesh%jjs(2,ms1)))**2)
3400  diff = sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - h_mesh%rr(:,h_mesh%jjs(1,ms2)))**2)
3401  IF (diff/ref .LT. 1.d-10) THEN ! 1 = 1
3402  w_cs = wws
3403  ELSE ! 1 = 2
3404  DO ls = 1, l_gs
3405  w_cs(1,ls)= wws(2,ls)
3406  w_cs(2,ls)= wws(1,ls)
3407  IF (n_ws1==3) w_cs(n_ws1,ls) = wws(n_ws1,ls)
3408  WRITE(*,*) ' Ouaps! oder of shape functions changed?'
3409  END DO
3410  END IF
3411 
3412  DO ls = 1, l_gs
3413  gauss2(1,ls) = sum(h_mesh%rr(1,h_mesh%jjs(:,ms2))*h_mesh%gauss%wws(:,ls))
3414  gauss2(2,ls) = sum(h_mesh%rr(2,h_mesh%jjs(:,ms2))*h_mesh%gauss%wws(:,ls))
3415  gauss1(1,ls) = sum(h_mesh%rr(1,h_mesh%jjs(:,ms1))*h_mesh%gauss%wws(:,ls))
3416  gauss1(2,ls) = sum(h_mesh%rr(2,h_mesh%jjs(:,ms1))*h_mesh%gauss%wws(:,ls))
3417  END DO
3418 
3419  DO ls2 = 1, l_gs
3420  ref = sqrt(1.d-8+sum(gauss2(:,ls2)**2))
3421  mark = .false.
3422  DO ls1 = 1, l_gs
3423  diff = sqrt(sum((gauss2(:,ls2)-gauss1(:,ls1))**2))
3424  IF (diff .LT. 1.d-10) THEN
3425  dw_cs(:,:,ls2,ms1) = h_mesh%gauss%dw_s(:,:,ls1,ms1)
3426  mark = .true.
3427  EXIT
3428  END IF
3429  END DO
3430  IF (.NOT.mark) WRITE(*,*) ' BUG '
3431  END DO
3432 
3433  END DO
3434 
3435  DO ms = 1, interface_h_mu%mes
3436 
3437  ms2 = interface_h_mu%mesh2(ms)
3438  ms1 = interface_h_mu%mesh1(ms)
3439  m2 = h_mesh%neighs(ms2)
3440  m1 = h_mesh%neighs(ms1)
3441  mu_h = sum(mu_h_field(h_mesh%jj(:,m1)))/h_mesh%gauss%n_w
3442  !JLG, FL, May, 28, 2009
3443  !hm1 = stab_colle_H_mu*(((mu_H+mu_H)/mu_H)/SUM(rjs(:,ms2)))
3444  hm1 = 1/sum(rjs(:,ms2))
3445  !JLG, FL, May, 28, 2009
3446  !====================================================================================
3447  !------------------------------------TERMES SUR LE BLOC H----------------------------
3448  !====================================================================================
3449 
3450  !-------------------------------hm1 (bi x ni) . (bj x nj)----------------------------
3451  !---------------------------------+ (mui bi.ni) (muj bj.nj)--------------------------
3452  !====================================================================================
3453  hsij = 0.d0
3454  DO ls = 1, l_gs
3455  !===Compute radius of Gauss point
3456  ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms2))* h_mesh%gauss%wws(:,ls))
3457  x = hm1*rjs(ls,ms2)*ray
3458 
3459  !June 14 2008, muhl
3460  muhl1 = sum(mu_h_field(h_mesh%jjs(:,ms1))*w_cs(:,ls))
3461  muhl2 = sum(mu_h_field(h_mesh%jjs(:,ms2))* wws(:,ls))
3462  !JLG, FL, May, 28, 2009, Normalization
3463  normt =stab_colle_h_mu
3464  norm = stab_div*sum(rjs(:,ms2))**(2*alpha)
3465  !norm = stab_div*SUM(rjs(:,ms2))**(2*alpha)/MAX(muhl1,muhl2)
3466  !norm = stab_div*SUM(rjs(:,ms2))**(2*alpha)/MAX(muhl1,muhl2)**2
3467  !norm = 1.d0/MAX(muhl1,muhl2)
3468  !norm = 1.d0/MAX(muhl1,muhl2)**2
3469  !JLG, FL, May, 28, 2009, Normalization
3470  !June 14 2008, muhl
3471 
3472  DO ci = 1, 2
3473  IF (ci==1) THEN
3474  normi = rnorms(:,ls,ms1)
3475  wwsi = w_cs(:,ls)
3476  n_wsi = n_ws1
3477  muhi = muhl1
3478  ELSE
3479  normi = rnorms(:,ls,ms2)
3480  wwsi = wws(:,ls)
3481  n_wsi = n_ws2
3482  muhi = muhl2
3483  END IF
3484  DO cj = 1, 2
3485  IF (cj==1) THEN
3486  normj = rnorms(:,ls,ms1)
3487  wwsj = w_cs(:,ls)
3488  n_wsj = n_ws1
3489  muhj = muhl1
3490  ELSE
3491  normj = rnorms(:,ls,ms2)
3492  wwsj = wws(:,ls)
3493  n_wsj = n_ws2
3494  muhj = muhl2
3495  END IF
3496 
3497  DO ni = 1, n_wsi
3498  DO nj = 1, n_wsj
3499  wwiwwj = x * wwsi(ni)*wwsj(nj)
3500  y = normt * wwiwwj
3501  ! June 14 2008, added z
3502  z = norm * muhi * muhj * wwiwwj
3503  ! June 14 2008, added z
3504  hsij(1,ni,nj,ci,cj) = hsij(1,ni,nj,ci,cj) + y*normi(2)*normj(2) &
3505  + z*normi(1)*normj(1)
3506  hsij(4,ni,nj,ci,cj) = hsij(4,ni,nj,ci,cj) - y*normj(1)*normi(2) &
3507  + z*normi(1)*normj(2)
3508  hsij(5,ni,nj,ci,cj) = hsij(5,ni,nj,ci,cj) + y*(normi(1)*normj(1) + normi(2)*normj(2))
3509  hsij(6,ni,nj,ci,cj) = hsij(6,ni,nj,ci,cj) + y*normi(1)*normj(1) &
3510  + z*normi(2)*normj(2)
3511  END DO
3512  END DO
3513  END DO
3514  END DO
3515  END DO
3516 
3517  mat_loc1 = 0.d0
3518  mat_loc2 = 0.d0
3519  DO ci = 1, 2
3520  DO ki = 1, 3
3521  DO ni = 1, n_ws1
3522  IF (ci==1) THEN
3523  i = interface_h_mu%jjs1(ni,ms)
3524  ELSE
3525  i = interface_h_mu%jjs2(ni,ms)
3526  END IF
3527  ib = la_h%loc_to_glob(ki,i)
3528  ix = ni + n_ws1*((ki-1) + 3*(ci-1))
3529  idxn(ix) = ib-1
3530 
3531  DO cj = 1, 2
3532  DO kj = 1, 3
3533  DO nj = 1, n_ws2
3534  IF (cj==1) THEN
3535  j = interface_h_mu%jjs1(nj,ms)
3536  ELSE
3537  j = interface_h_mu%jjs2(nj,ms)
3538  END IF
3539  jb = la_h%loc_to_glob(kj,j)
3540  jx = nj + n_ws2*((kj-1) + 3*(cj-1))
3541  jdxn(jx) = jb-1
3542  IF ((ki == 1) .AND. (kj == 1)) THEN
3543  mat_loc1(ix,jx) = hsij(1,ni,nj,ci,cj)
3544  mat_loc2(ix,jx) = hsij(1,ni,nj,ci,cj)
3545  ELSEIF ((ki == 1) .AND. (kj == 3)) THEN
3546  mat_loc1(ix,jx) = hsij(4,ni,nj,ci,cj)
3547  mat_loc2(ix,jx) = hsij(4,ni,nj,ci,cj)
3548  ELSEIF ((ki == 3) .AND. (kj == 1)) THEN
3549  mat_loc1(ix,jx) = hsij(4,nj,ni,cj,ci)
3550  mat_loc2(ix,jx) = hsij(4,nj,ni,cj,ci)
3551  ELSEIF ((ki == 2) .AND. (kj == 2)) THEN
3552  mat_loc1(ix,jx) = hsij(5,ni,nj,ci,cj)
3553  mat_loc2(ix,jx) = hsij(5,ni,nj,ci,cj)
3554  ELSEIF ((ki == 3) .AND. (kj == 3)) THEN
3555  mat_loc1(ix,jx) = hsij(6,ni,nj,ci,cj)
3556  mat_loc2(ix,jx) = hsij(6,ni,nj,ci,cj)
3557  ENDIF
3558  END DO
3559  END DO
3560  END DO
3561  END DO
3562  END DO
3563  END DO
3564 
3565  CALL matsetvalues(h_p_phi_mat1, 6*n_ws1, idxn(1:6*n_ws1), 6*n_ws2, jdxn(1:6*n_ws2), &
3566  mat_loc1(1:6*n_ws1,1:6*n_ws2), add_values, ierr)
3567  CALL matsetvalues(h_p_phi_mat2, 6*n_ws1, idxn(1:6*n_ws1), 6*n_ws2, jdxn(1:6*n_ws2), &
3568  mat_loc2(1:6*n_ws1,1:6*n_ws2), add_values, ierr)
3569 
3570  !====================================================================================
3571  !------------------------(1/sigma) (Rot bj) . (bi x ni)------------------------------
3572  !====================================================================================
3573 
3574  !terme sans derivee
3575  hsij = 0.d0
3576  gsij = 0.d0
3577  DO ls = 1, h_mesh%gauss%l_Gs
3578 
3579  !===Compute radius of Gauss point
3580  ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms2))* h_mesh%gauss%wws(:,ls))
3581  x = rjs(ls,ms2)*ray
3582 ! TEST DEBUG
3583 
3584 ! TEST DEBUG
3585  DO ci = 1, 2
3586  IF (ci==1) THEN
3587  normi = rnorms(:,ls,ms1)
3588  wwsi = w_cs(:,ls)
3589  n_wsi = n_ws1
3590  sigmai = sigma(m1)
3591 ! TEST DEBUG
3592 
3593 ! TEST DEBUG
3594  ELSE
3595  normi = rnorms(:,ls,ms2)
3596  wwsi = wws(:,ls)
3597  n_wsi = n_ws2
3598  sigmai = sigma(m2)
3599 ! TEST DEBUG
3600 
3601 ! TEST DEBUG
3602  END IF
3603  DO cj = 1, 2
3604  IF (cj==1) THEN
3605  normj = rnorms(:,ls,ms1)
3606  wwsj = w_cs(:,ls)
3607  n_wsj = n_ws1
3608  sigmaj = sigma(m1)
3609 ! TEST DEBUG
3610 
3611 ! TEST DEBUG
3612  ELSE
3613  normj = rnorms(:,ls,ms2)
3614  wwsj = wws(:,ls)
3615  n_wsj = n_ws2
3616  sigmaj = sigma(m2)
3617 ! TEST DEBUG
3618 
3619 ! TEST DEBUG
3620  END IF
3621 
3622  DO ni = 1,n_wsi !
3623  DO nj = 1, n_wsj!
3624  y = x*wwsi(ni)*wwsj(nj)/(2*sigmaj)
3625  hsij(2,ni,nj,ci,cj) = hsij(2,ni,nj,ci,cj) + y * (-mode/ray)*normi(1)
3626  hsij(3,ni,nj,ci,cj) = hsij(3,ni,nj,ci,cj) + y * mode/ray *normi(1)
3627  hsij(5,ni,nj,ci,cj) = hsij(5,ni,nj,ci,cj) + y * (-1/ray) *normi(1)
3628  hsij(8,ni,nj,ci,cj) = hsij(8,ni,nj,ci,cj) + y * (-mode/ray)*normi(2)
3629  hsij(9,ni,nj,ci,cj) = hsij(9,ni,nj,ci,cj) + y * mode/ray *normi(2)
3630  y = x*wwsi(ni)*wwsj(nj)/(2*sigmai)
3631  gsij(2,ni,nj,ci,cj) = gsij(2,ni,nj,ci,cj) + y * (-mode/ray)*normj(1)
3632  gsij(3,ni,nj,ci,cj) = gsij(3,ni,nj,ci,cj) + y * ( mode/ray)*normj(1)
3633  gsij(5,ni,nj,ci,cj) = gsij(5,ni,nj,ci,cj) + y * (-1/ray) *normj(1)
3634  gsij(8,ni,nj,ci,cj) = gsij(8,ni,nj,ci,cj) + y * (-mode/ray)*normj(2)
3635  gsij(9,ni,nj,ci,cj) = gsij(9,ni,nj,ci,cj) + y * mode/ray *normj(2)
3636  ENDDO
3637  ENDDO
3638  ENDDO
3639  END DO
3640  END DO
3641 
3642  !June 14 2008
3643  gsij = c_sym*gsij
3644  !June 14 2008
3645 
3646  mat_loc1 = 0.d0
3647  mat_loc2 = 0.d0
3648 
3649  DO ci = 1, 2
3650  DO ki = 1, 3
3651  DO ni = 1, n_wsi
3652  IF (ci==1) THEN
3653  i = interface_h_mu%jjs1(ni,ms)
3654  ELSE
3655  i = interface_h_mu%jjs2(ni,ms)
3656  END IF
3657  ib = la_h%loc_to_glob(ki,i)
3658  ix = ni + n_wsi*((ki-1) + 3*(ci-1))
3659  idxn(ix) = ib - 1
3660  DO cj = 1, 2
3661  DO kj = 1, 3
3662  DO nj = 1, n_wsj
3663  IF (cj==1) THEN
3664  j = interface_h_mu%jjs1(nj,ms)
3665  ELSE
3666  j = interface_h_mu%jjs2(nj,ms)
3667  END IF
3668  jb = la_h%loc_to_glob(kj,j)
3669  jx = nj + n_wsj*((kj-1) + 3*(cj-1))
3670  jdxn(jx) = jb - 1
3671  IF ((ki == 2) .AND. (kj == 1)) THEN
3672  mat_loc1(ix,jx) = hsij(2,ni,nj,ci,cj)
3673  mat_loc2(ix,jx) = hsij(3,ni,nj,ci,cj)
3674  ELSE IF((ki == 1) .AND. (kj == 2)) THEN
3675  mat_loc1(ix,jx) = gsij(2,ni,nj,ci,cj)
3676  mat_loc2(ix,jx) = gsij(3,ni,nj,ci,cj)
3677  ELSEIF ((ki == 2) .AND. (kj == 2)) THEN
3678  mat_loc1(ix,jx) = hsij(5,ni,nj,ci,cj)+gsij(5,ni,nj,ci,cj)
3679  mat_loc2(ix,jx) = hsij(5,ni,nj,ci,cj)+gsij(5,ni,nj,ci,cj)
3680  ELSEIF ((ki == 2) .AND. (kj == 3)) THEN
3681  mat_loc1(ix,jx) = hsij(8,ni,nj,ci,cj)
3682  mat_loc2(ix,jx) = hsij(9,ni,nj,ci,cj)
3683  ELSEIF ((ki == 3) .AND. (kj == 2)) THEN
3684  mat_loc1(ix,jx) = gsij(8,ni,nj,ci,cj)
3685  mat_loc2(ix,jx) = gsij(9,ni,nj,ci,cj)
3686  ENDIF
3687  END DO
3688  END DO
3689  END DO
3690  END DO
3691  END DO
3692  END DO
3693 
3694  CALL matsetvalues(h_p_phi_mat1, 6*n_ws1, idxn(1:6*n_ws1), 6*n_ws2, jdxn(1:6*n_ws2), &
3695  mat_loc1(1:6*n_ws1,1:6*n_ws2), add_values, ierr)
3696  CALL matsetvalues(h_p_phi_mat2, 6*n_ws1, idxn(1:6*n_ws1), 6*n_ws2, jdxn(1:6*n_ws2), &
3697  mat_loc2(1:6*n_ws1,1:6*n_ws2), add_values, ierr)
3698 
3699  !terme avec derivees
3700  hsij = 0.d0
3701  gsij = 0.d0
3702  DO ls = 1, h_mesh%gauss%l_Gs
3703 
3704  !===Compute radius of Gauss point
3705  ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms2))* h_mesh%gauss%wws(:,ls))
3706  x = rjs(ls,ms2)*ray
3707 ! TEST DEBUG
3708 
3709 ! TEST DEBUG
3710 
3711  DO ci = 1, 2
3712  IF (ci==1) THEN
3713  normi = rnorms(:,ls,ms1)
3714  wwsi = w_cs(:,ls)
3715  dwsi = dw_cs(:,:,ls,ms1)
3716  n_wsi = n_ws1
3717  n_wi = n_w1
3718  sigmai = sigma(m1)
3719 ! TEST DEBUG
3720 
3721 ! TEST DEBUG
3722  ELSE
3723  normi = rnorms(:,ls,ms2)
3724  wwsi = wws(:,ls)
3725  dwsi = dw_s(:,:,ls,ms2)
3726  n_wsi = n_ws2
3727  n_wi = n_w2
3728  sigmai = sigma(m2)
3729 ! TEST DEBUG
3730 
3731 ! TEST DEBUG
3732  END IF
3733  DO cj = 1, 2
3734  IF (cj==1) THEN
3735  normj = rnorms(:,ls,ms1)
3736  wwsj = w_cs(:,ls)
3737  dwsj = dw_cs(:,:,ls,ms1)
3738  n_wsj = n_ws1
3739  n_wj = n_w1
3740  sigmaj = sigma(m1)
3741 ! TEST DEBUG
3742 
3743 ! TEST DEBUG
3744  ELSE
3745  normj = rnorms(:,ls,ms2)
3746  wwsj = wws(:,ls)
3747  dwsj = dw_s(:,:,ls,ms2)
3748  n_wsj = n_ws2
3749  n_wj = n_w2
3750  sigmaj = sigma(m2)
3751 ! TEST DEBUG
3752 
3753 ! TEST DEBUG
3754  END IF
3755 
3756  !termes avec derivees
3757  DO ni = 1,n_wsi
3758  DO nj = 1, n_wj
3759  y = x*wwsi(ni)/(2*sigmaj)
3760  hsij(1,ni,nj,ci,cj) = hsij(1,ni,nj,ci,cj) + y*(-dwsj(2,nj))*normi(2)
3761  hsij(4,ni,nj,ci,cj) = hsij(4,ni,nj,ci,cj) + y* dwsj(1,nj) *normi(2)
3762  hsij(5,ni,nj,ci,cj) = hsij(5,ni,nj,ci,cj) + y*(-dwsj(2,nj) *normi(2)-dwsj(1,nj)*normi(1))
3763  hsij(6,ni,nj,ci,cj) = hsij(6,ni,nj,ci,cj) + y*(-dwsj(1,nj))*normi(1)
3764  hsij(7,ni,nj,ci,cj) = hsij(7,ni,nj,ci,cj) + y* dwsj(2,nj) *normi(1)
3765  ENDDO
3766  END DO
3767  DO ni = 1,n_wi
3768  DO nj = 1, n_wsj
3769  y = x*wwsj(nj)/(2*sigmai)
3770  gsij(1,ni,nj,ci,cj) = gsij(1,ni,nj,ci,cj) + y*(-dwsi(2,ni))*normj(2)
3771  gsij(4,ni,nj,ci,cj) = gsij(4,ni,nj,ci,cj) + y* dwsi(2,ni) *normj(1)
3772  gsij(5,ni,nj,ci,cj) = gsij(5,ni,nj,ci,cj) + y*(-dwsi(2,ni) *normj(2)-dwsi(1,ni)*normj(1))
3773  gsij(6,ni,nj,ci,cj) = gsij(6,ni,nj,ci,cj) + y*(-dwsi(1,ni))*normj(1)
3774  gsij(7,ni,nj,ci,cj) = gsij(7,ni,nj,ci,cj) + y* dwsi(1,ni) *normj(2)
3775  ENDDO
3776  END DO
3777 
3778  ENDDO
3779  ENDDO
3780  ENDDO
3781 
3782  !June 14 2008
3783  gsij = c_sym*gsij
3784  !June 14 2008
3785 
3786  mat_loc1 = 0.d0
3787  mat_loc2 = 0.d0
3788  DO ci = 1, 2
3789  DO ki = 1, 3
3790  DO ni = 1, n_wsi
3791  IF (ci==1) THEN
3792  i = interface_h_mu%jjs1(ni,ms)
3793  ELSE
3794  i = interface_h_mu%jjs2(ni,ms)
3795  END IF
3796  ib = la_h%loc_to_glob(ki,i)
3797  ix = ni + n_wsi*((ki-1) + 3*(ci-1))
3798  idxn(ix) = ib - 1
3799  DO cj = 1, 2
3800  DO kj = 1, 3
3801  DO nj = 1, n_wj
3802  IF (cj==1) THEN
3803  j = h_mesh%jj(nj,m1)
3804  ELSE
3805  j = h_mesh%jj(nj,m2)
3806  END IF
3807  jb = la_h%loc_to_glob(kj,j)
3808  jx = nj + n_wj*((kj-1) + 3*(cj-1))
3809  jdxn(jx) = jb - 1
3810  IF ((ki == 1) .AND. (kj == 1)) THEN
3811  mat_loc1(ix,jx) = hsij(1,ni,nj,ci,cj)
3812  mat_loc2(ix,jx) = hsij(1,ni,nj,ci,cj)
3813  ELSEIF ((ki == 1) .AND. (kj == 3)) THEN
3814  mat_loc1(ix,jx) = hsij(4,ni,nj,ci,cj)
3815  mat_loc2(ix,jx) = hsij(4,ni,nj,ci,cj)
3816  ELSEIF ((ki == 3) .AND. (kj == 1)) THEN
3817  mat_loc1(ix,jx) = hsij(7,ni,nj,ci,cj)
3818  mat_loc2(ix,jx) = hsij(7,ni,nj,ci,cj)
3819  ELSEIF ((ki == 2) .AND. (kj == 2)) THEN
3820  mat_loc1(ix,jx) = hsij(5,ni,nj,ci,cj)
3821  mat_loc2(ix,jx) = hsij(5,ni,nj,ci,cj)
3822  ELSEIF ((ki == 3) .AND. (kj == 3)) THEN
3823  mat_loc1(ix,jx) = hsij(6,ni,nj,ci,cj)
3824  mat_loc2(ix,jx) = hsij(6,ni,nj,ci,cj)
3825  ENDIF
3826 
3827  END DO
3828  END DO
3829  END DO
3830  END DO
3831  END DO
3832  END DO
3833 
3834  CALL matsetvalues(h_p_phi_mat1, 6*n_ws1, idxn(1:6*n_ws1), 6*n_w2, jdxn(1:6*n_w2), &
3835  mat_loc1(1:6*n_ws1,1:6*n_w2), add_values, ierr)
3836  CALL matsetvalues(h_p_phi_mat2, 6*n_ws1, idxn(1:6*n_ws1), 6*n_w2, jdxn(1:6*n_w2), &
3837  mat_loc2(1:6*n_ws1,1:6*n_w2), add_values, ierr)
3838 
3839  mat_loc1 = 0.d0
3840  mat_loc2 = 0.d0
3841  DO ci = 1, 2
3842  DO ki = 1, 3
3843  DO ni = 1, n_wi
3844  IF (ci==1) THEN
3845  i = h_mesh%jj(ni,m1)
3846  ELSE
3847  i = h_mesh%jj(ni,m2)
3848  END IF
3849  ib = la_h%loc_to_glob(ki,i)
3850  ix = ni + n_wi*((ki-1) + 3*(ci-1))
3851  idxn(ix) = ib - 1
3852  DO cj = 1, 2
3853  DO kj = 1, 3
3854  DO nj = 1, n_wsj
3855  IF (cj==1) THEN
3856  j = interface_h_mu%jjs1(nj,ms)
3857  ELSE
3858  j = interface_h_mu%jjs2(nj,ms)
3859  END IF
3860  jb = la_h%loc_to_glob(kj,j)
3861  jx = nj + n_wsj*((kj-1) + 3*(cj-1))
3862  jdxn(jx) = jb-1
3863  IF ((ki == 1) .AND. (kj == 1)) THEN
3864  mat_loc1(ix,jx) = gsij(1,ni,nj,ci,cj)
3865  mat_loc2(ix,jx) = gsij(1,ni,nj,ci,cj)
3866  ELSEIF ((ki == 1) .AND. (kj == 3)) THEN
3867  mat_loc1(ix,jx) = gsij(4,ni,nj,ci,cj)
3868  mat_loc2(ix,jx) = gsij(4,ni,nj,ci,cj)
3869  ELSEIF ((ki == 3) .AND. (kj == 1)) THEN
3870  mat_loc1(ix,jx) = gsij(7,ni,nj,ci,cj)
3871  mat_loc2(ix,jx) = gsij(7,ni,nj,ci,cj)
3872  ELSEIF ((ki == 2) .AND. (kj == 2)) THEN
3873  mat_loc1(ix,jx) = gsij(5,ni,nj,ci,cj)
3874  mat_loc2(ix,jx) = gsij(5,ni,nj,ci,cj)
3875  ELSEIF ((ki == 3) .AND. (kj == 3)) THEN
3876  mat_loc1(ix,jx) = gsij(6,ni,nj,ci,cj)
3877  mat_loc2(ix,jx) = gsij(6,ni,nj,ci,cj)
3878  ENDIF
3879  END DO
3880  END DO
3881  END DO
3882  END DO
3883  END DO
3884  END DO
3885 
3886  CALL matsetvalues(h_p_phi_mat1, 6*n_w1, idxn(1:6*n_w1), 6*n_ws2, jdxn(1:6*n_ws2), &
3887  mat_loc1(1:6*n_w1,1:6*n_ws2), add_values, ierr)
3888  CALL matsetvalues(h_p_phi_mat2, 6*n_w1, idxn(1:6*n_w1), 6*n_ws2, jdxn(1:6*n_ws2), &
3889  mat_loc2(1:6*n_w1,1:6*n_ws2), add_values, ierr)
3890 
3891  END DO
3892 
3893  CALL matassemblybegin(h_p_phi_mat1,mat_final_assembly,ierr)
3894  CALL matassemblyend(h_p_phi_mat1,mat_final_assembly,ierr)
3895  CALL matassemblybegin(h_p_phi_mat2,mat_final_assembly,ierr)
3896  CALL matassemblyend(h_p_phi_mat2,mat_final_assembly,ierr)
3897 
3898 !!$ DEALLOCATE(mat_loc1, mat_loc2, idxn, jdxn)
3899 !!$ DEALLOCATE(w_cs, dw_cs, gauss1, gauss2, dwsi, dwsj)
3900 
3901  END SUBROUTINE mat_maxwell_mu
3902 
3903  SUBROUTINE courant_mu(H_mesh,interface_H_mu,sigma,mu_H_field,time,mode,nl, &
3904  la_h, vb_1, vb_2, h_ext)
3905  !forcage faisant intervenir J, volumique et interface pour H et phi
3906  !pour le probleme en entier
3907 
3908  USE def_type_mesh
3909  USE gauss_points
3910  USE boundary
3911 
3912  IMPLICIT NONE
3913  TYPE(mesh_type), INTENT(IN) :: h_mesh
3914  TYPE(interface_type), INTENT(IN) :: interface_h_mu
3915  REAL(KIND=8), INTENT(IN) :: time
3916  REAL(KIND=8), DIMENSION(H_mesh%me), INTENT(IN) :: sigma
3917  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: mu_h_field
3918  INTEGER, INTENT(IN) :: mode
3919  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: nl
3920  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: h_ext
3921  REAL(KIND=8), DIMENSION(H_mesh%np,6) :: src_h
3922 
3923  REAL(KIND=8), DIMENSION(H_mesh%gauss%n_ws,H_mesh%gauss%l_Gs) :: w_cs
3924  REAL(KIND=8), DIMENSION(2) :: normi, gaussp1, gaussp2
3925  REAL(KIND=8), DIMENSION(H_mesh%gauss%n_ws) :: wwsi
3926  REAL(KIND=8) :: x, ray
3927  INTEGER :: i, ni, ms, k, ls, n_ws1, n_ws2, ms1, ms2, n_w1, n_w2, m1, m2, ci, n_wsi
3928  INTEGER :: mesh_id1, mesh_id2
3929  REAL(KIND=8), DIMENSION(6) :: jsolh_anal, test, h_ext_l
3930  REAL(KIND=8) :: muhl1, muhl2, ref, diff
3931  !April 17th, 2008, JLG
3932  REAL(KIND=8) :: one
3933  DATA one/1.d0/
3934  !April 17th, 2008, JLG
3935 !$$ FL+CN 22/03/2013
3936 !!$ INTEGER, DIMENSION(:), ALLOCATABLE :: idxn
3937  INTEGER, DIMENSION(H_mesh%np) :: idxn
3938 !$$ FL+CN 22/03/2013
3939  TYPE(petsc_csr_la) :: la_h
3940 #include "petsc/finclude/petsc.h"
3941  petscerrorcode :: ierr
3942  vec :: vb_1, vb_2
3943 
3944 
3945  !**********************************************************************************
3946  !--------------------TERMS ON SIGMA_MU-------------------------------
3947  !**********************************************************************************
3948  src_h = 0.d0
3949  !WRITE(*,*) 'Assembling rhs interface_H_mu'
3950  CALL gauss(h_mesh)
3951  n_ws1 = h_mesh%gauss%n_ws
3952  n_ws2 = h_mesh%gauss%n_ws
3953  n_w1 = h_mesh%gauss%n_w
3954  n_w2 = h_mesh%gauss%n_w
3955 
3956  DO ms = 1, interface_h_mu%mes
3957  ms1 = interface_h_mu%mesh1(ms)
3958  ms2 = interface_h_mu%mesh2(ms)
3959  m1 = h_mesh%neighs(ms1)
3960  m2 = h_mesh%neighs(ms2)
3961 
3962  ref = 1.d-8+sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - h_mesh%rr(:,h_mesh%jjs(2,ms1)))**2)
3963  diff = sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - h_mesh%rr(:,h_mesh%jjs(1,ms2)))**2)
3964  IF (diff/ref .LT. 1.d-10) THEN ! 1 = 1
3965  w_cs = wws
3966  ELSE ! 1 = 2
3967  DO ls = 1, l_gs
3968  w_cs(1,ls)= wws(2,ls)
3969  w_cs(2,ls)= wws(1,ls)
3970  IF (n_ws1==3) w_cs(n_ws1,ls) = wws(n_ws1,ls)
3971  WRITE(*,*) ' Ouaps! oder of shape functions changed?'
3972  END DO
3973  END IF
3974  END DO
3975 
3976  DO ms = 1, interface_h_mu%mes
3977  ms2 = interface_h_mu%mesh2(ms)
3978  ms1 = interface_h_mu%mesh1(ms)
3979  m2 = h_mesh%neighs(ms2)
3980  m1 = h_mesh%neighs(ms1)
3981  mesh_id1 = h_mesh%i_d(m1)
3982  mesh_id2 = h_mesh%i_d(m2)
3983  DO ls = 1, l_gs
3984  !===Compute radius of Gauss point
3985  ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms2))* h_mesh%gauss%wws(:,ls))
3986 
3987 
3988  ! Cote 1
3989  DO k=1, 6
3990  h_ext_l(k) = sum(h_ext(h_mesh%jjs(:,ms1),k)*h_mesh%gauss%wws(:,ls))
3991  END DO
3992  muhl1=sum(mu_h_field(h_mesh%jjs(:,ms1))*w_cs(:,ls))
3993  gaussp1(1) = sum(h_mesh%rr(1,h_mesh%jjs(:,ms1))*w_cs(:,ls))
3994  gaussp1(2) = sum(h_mesh%rr(2,h_mesh%jjs(:,ms1))*w_cs(:,ls))
3995  DO k=1, 6
3996  jsolh_anal(k) = jexact_gauss(k, gaussp1, mode, one ,sigma(m1), &
3997  muhl1, time, mesh_id1, h_ext_l)/sigma(m1) &
3998  + muhl1 * sum(nl(h_mesh%jjs(1:n_ws1,ms1),k)*w_cs(1:n_ws1,ls))
3999 ! TEST DEBUG
4000 
4001 ! TEST DEBUG
4002  ENDDO
4003 
4004  ! Cote 2
4005  DO k=1, 6
4006  h_ext_l(k) = sum(h_ext(h_mesh%jjs(:,ms2),k)*h_mesh%gauss%wws(:,ls))
4007  END DO
4008  muhl2=sum(mu_h_field(h_mesh%jjs(:,ms2))*wws(:,ls))
4009  gaussp2(1) = sum(h_mesh%rr(1,h_mesh%jjs(:,ms2))*wws(:,ls))
4010  gaussp2(2) = sum(h_mesh%rr(2,h_mesh%jjs(:,ms2))*wws(:,ls))
4011  IF (maxval(abs(gaussp1-gaussp2)) > 1.d-11) THEN
4012  WRITE(*,*) ' BUG courant_mu '
4013  stop
4014  END IF
4015  DO k=1, 6
4016  test(k) = jexact_gauss(k, gaussp2, mode, one ,sigma(m2), &
4017  muhl2, time, mesh_id2, h_ext_l)/sigma(m2) &
4018  + muhl2 * sum(nl(h_mesh%jjs(1:n_ws2,ms2),k)*wws(1:n_ws2,ls))
4019  !JsolH_anal(k) = JsolH_anal(k) + Jexact_gauss(k, gaussp2, mode, one,sigma(m2), muhl2, time)/sigma(m2) &
4020  ! + muhl2 * SUM(NL(H_mesh%jjs(1:n_ws2,ms2),k)*wws(1:n_ws2,ls))
4021  !WRITE(*,*) ABS(test(k) - JsolH_anal(k)), test(k), JsolH_anal(k)
4022  jsolh_anal(k) = jsolh_anal(k) + test(k)
4023 ! TEST DEBUG
4024 
4025 ! TEST DEBUG
4026  ENDDO
4027  ! Division by 2 to get the mean is in definition of x below.
4028 
4029  !---------forcage pour H
4030  DO ci = 1, 2
4031  IF (ci==1) THEN
4032  normi = rnorms(:,ls,ms1)
4033  wwsi = w_cs(:,ls)
4034  n_wsi = n_ws1
4035  ELSE
4036  normi = rnorms(:,ls,ms2)
4037  wwsi = wws(:,ls)
4038  n_wsi = n_ws2
4039  END IF
4040  DO ni = 1, n_wsi
4041  IF (ci==1) THEN
4042  i = interface_h_mu%jjs1(ni,ms)
4043  ELSE
4044  i = interface_h_mu%jjs2(ni,ms)
4045  END IF
4046  x = rjs(ls,ms2)*ray*wwsi(ni)/2
4047  src_h(i,1) = src_h(i,1)+x*(-jsolh_anal(3)*normi(2))
4048  src_h(i,2) = src_h(i,2)+x*(-jsolh_anal(4)*normi(2))
4049  src_h(i,3) = src_h(i,3)+x*(jsolh_anal(1)*normi(2)-jsolh_anal(5)*normi(1))
4050  src_h(i,4) = src_h(i,4)+x*(jsolh_anal(2)*normi(2)-jsolh_anal(6)*normi(1))
4051  src_h(i,5) = src_h(i,5)+x*(jsolh_anal(3)*normi(1))
4052  src_h(i,6) = src_h(i,6)+x*(jsolh_anal(4)*normi(1))
4053  END DO
4054  ENDDO
4055  END DO
4056  END DO
4057 
4058  IF (h_mesh%np /= 0) THEN
4059 !!$ ALLOCATE(idxn(H_mesh%np))
4060  idxn = la_h%loc_to_glob(1,:)-1
4061  CALL vecsetvalues(vb_1, h_mesh%np, idxn, src_h(:,1), add_values, ierr)
4062  CALL vecsetvalues(vb_2, h_mesh%np, idxn, src_h(:,2), add_values, ierr)
4063  idxn = la_h%loc_to_glob(2,:)-1
4064  CALL vecsetvalues(vb_1, h_mesh%np, idxn, src_h(:,4), add_values, ierr)
4065  CALL vecsetvalues(vb_2, h_mesh%np, idxn, src_h(:,3), add_values, ierr)
4066  idxn = la_h%loc_to_glob(3,:)-1
4067  CALL vecsetvalues(vb_1, h_mesh%np, idxn, src_h(:,5), add_values, ierr)
4068  CALL vecsetvalues(vb_2, h_mesh%np, idxn, src_h(:,6), add_values, ierr)
4069 !!$ DEALLOCATE(idxn)
4070  END IF
4071  CALL vecassemblybegin(vb_1,ierr)
4072  CALL vecassemblyend(vb_1,ierr)
4073  CALL vecassemblybegin(vb_2,ierr)
4074  CALL vecassemblyend(vb_2,ierr)
4075 
4076  END SUBROUTINE courant_mu
4077 
4078  SUBROUTINE rhs_dirichlet(H_mesh,Dirichlet_bdy_H_sides,sigma,&
4079  mu_h_field,time,mode,nl,stab, la_h, vb_1, vb_2, h_ext, j_over_sigma, sigma_curl_bdy)
4080  !forcage faisant intervenir J, volumique et surfacique
4081  !pour le probleme en entier
4082  USE def_type_mesh
4083  USE boundary
4084  USE input_data
4085 
4086  IMPLICIT NONE
4087  TYPE(mesh_type), INTENT(IN) :: h_mesh
4088  INTEGER, DIMENSION(:), INTENT(IN) :: dirichlet_bdy_h_sides
4089  REAL(KIND=8), INTENT(IN) :: time
4090  REAL(KIND=8), DIMENSION(H_mesh%me), INTENT(IN) :: sigma
4091  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: mu_h_field
4092  INTEGER, INTENT(IN) :: mode
4093  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: nl
4094  REAL(KIND=8), DIMENSION(3), INTENT(IN) :: stab
4095  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: h_ext
4096  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: j_over_sigma !Used only if sigma variable in fluid
4097  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: sigma_curl_bdy !Used only if sigma variable in fluid
4098  REAL(KIND=8), DIMENSION(H_mesh%np,6) :: src_h
4099  REAL(KIND=8), DIMENSION(2) :: gaussp1
4100  REAL(KIND=8) :: x, ray, stab_colle_h_mu
4101  INTEGER :: i, ni, ms, k, ls, m1, count
4102  INTEGER :: mesh_id1
4103  REAL(KIND=8), DIMENSION(6) :: jsolh_anal, h_ext_l
4104  REAL(KIND=8) :: muhl1, hm1
4105  REAL(KIND=8), DIMENSION(6,H_mesh%gauss%l_Gs) :: hloc, hlocxn
4106  REAL(KIND=8), DIMENSION(2,H_mesh%gauss%l_Gs) :: rloc
4107  REAL(KIND=8), DIMENSION(1) :: muloc
4108  INTEGER :: index
4109  !April 17th, 2008, JLG
4110  REAL(KIND=8) :: one
4111  DATA one/1.d0/
4112  !April 17th, 2008, JLG
4113 !!$ FL+CN 22/03/2013
4114 !!$ INTEGER, DIMENSION(:), ALLOCATABLE :: idxn
4115  INTEGER, DIMENSION(H_mesh%np) :: idxn
4116 !!$ FL+CN 22/03/2013
4117  TYPE(petsc_csr_la) :: la_h
4118 #include "petsc/finclude/petsc.h"
4119  petscerrorcode :: ierr
4120  vec :: vb_1, vb_2
4121 
4122  !IF (SIZE(Dirichlet_bdy_H_sides)==0) THEN
4123  ! RETURN
4124  !END IF
4125 
4126  src_h = 0.d0
4127 
4128  !IF (SIZE(Dirichlet_bdy_H_sides)==0) THEN
4129  ! IF (ASSOCIATED(nl)) DEALLOCATE(nl)
4130  ! RETURN
4131  !END IF
4132 
4133  stab_colle_h_mu = stab(3)
4134  index = 0
4135 
4136  DO count = 1, SIZE(dirichlet_bdy_h_sides)
4137  ms = dirichlet_bdy_h_sides(count)
4138  hm1 = stab_colle_h_mu/sum(h_mesh%gauss%rjs(:,ms))
4139  m1 = h_mesh%neighs(ms)
4140  mesh_id1 = h_mesh%i_d(m1)
4141  muloc(1) = mu_h_field(h_mesh%jj(1,m1))
4142 
4143  DO ls = 1, h_mesh%gauss%l_Gs
4144  rloc(1,ls) = sum(h_mesh%rr(1,h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
4145  rloc(2,ls) = sum(h_mesh%rr(2,h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
4146  END DO
4147 
4148  DO k = 1, 6
4149  hloc(k,:) = hexact(h_mesh, k, rloc, mode, muloc, time)
4150  END DO
4151 
4152  hlocxn(1,:) = hloc(3,:)*h_mesh%gauss%rnorms(2,:,ms)
4153  hlocxn(2,:) = hloc(4,:)*h_mesh%gauss%rnorms(2,:,ms)
4154  hlocxn(3,:) = hloc(5,:)*h_mesh%gauss%rnorms(1,:,ms)-hloc(1,:)*h_mesh%gauss%rnorms(2,:,ms)
4155  hlocxn(4,:) = hloc(6,:)*h_mesh%gauss%rnorms(1,:,ms)-hloc(2,:)*h_mesh%gauss%rnorms(2,:,ms)
4156  hlocxn(5,:) = -hloc(3,:)*h_mesh%gauss%rnorms(1,:,ms)
4157  hlocxn(6,:) = -hloc(4,:)*h_mesh%gauss%rnorms(1,:,ms)
4158 
4159  DO ls = 1, h_mesh%gauss%l_Gs
4160  index = index + 1
4161  !===Compute radius of Gauss point
4162  ray = rloc(1,ls) !SUM(H_mesh%rr(1,H_mesh%jjs(:,ms))* H_mesh%gauss%wws(:,ls))
4163  DO k = 1, 6
4164  h_ext_l(k) = sum(h_ext(h_mesh%jjs(:,ms),k)*h_mesh%gauss%wws(:,ls))
4165  END DO
4166 
4167  ! Cote 1
4168  muhl1=sum(mu_h_field(h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
4169  gaussp1(1) = rloc(1,ls) !SUM(H_mesh%rr(1,H_mesh%jjs(:,ms))*H_mesh%gauss%wws(:,ls))
4170  gaussp1(2) = rloc(2,ls) !SUM(H_mesh%rr(2,H_mesh%jjs(:,ms))*H_mesh%gauss%wws(:,ls))
4171 !!$ DO k=1, 6
4172 !!$ JsolH_anal(k) = Jexact_gauss(k, gaussp1, mode, one ,sigma(m1), muhl1, &
4173 !!$ time, mesh_id1, H_ext_l)/sigma(m1) &
4174 !!$ + muhl1 * SUM(NL(H_mesh%jjs(:,ms),k)*H_mesh%gauss%wws(:,ls)) &
4175 !!$ + hm1*Hlocxn(k,ls)
4176 !!$ ENDDO
4177 ! TEST DEBUG
4178  IF (inputs%if_level_set.AND.inputs%variation_sigma_fluid) THEN
4179  DO k = 1, 6
4180  jsolh_anal(k) = sum(j_over_sigma(h_mesh%jjs(:,ms),k)*h_mesh%gauss%wws(:,ls)) &
4181  + sigma_curl_bdy(index,k) &
4182  + sum(nl(h_mesh%jjs(:,ms),k)*h_mesh%gauss%wws(:,ls)) &
4183  + hm1*hlocxn(k,ls)
4184  END DO
4185  ELSE
4186  DO k = 1, 6
4187  jsolh_anal(k) = jexact_gauss(k, gaussp1, mode, one ,sigma(m1), muhl1, &
4188  time, mesh_id1, h_ext_l)/sigma(m1) &
4189  + sum(nl(h_mesh%jjs(:,ms),k)*h_mesh%gauss%wws(:,ls)) &
4190  + hm1*hlocxn(k,ls)
4191  END DO
4192  END IF
4193 ! TEST DEBUG
4194 
4195  DO ni = 1, h_mesh%gauss%n_ws
4196  i = h_mesh%jjs(ni,ms)
4197  x = h_mesh%gauss%rjs(ls,ms)*ray*h_mesh%gauss%wws(ni,ls)
4198 
4199  src_h(i,1) = src_h(i,1)+x*(-jsolh_anal(3)*h_mesh%gauss%rnorms(2,ls,ms))
4200  src_h(i,2) = src_h(i,2)+x*(-jsolh_anal(4)*h_mesh%gauss%rnorms(2,ls,ms))
4201  src_h(i,3) = src_h(i,3)+x*(jsolh_anal(1)*h_mesh%gauss%rnorms(2,ls,ms)&
4202  -jsolh_anal(5)*h_mesh%gauss%rnorms(1,ls,ms))
4203  src_h(i,4) = src_h(i,4)+x*(jsolh_anal(2)*h_mesh%gauss%rnorms(2,ls,ms)&
4204  -jsolh_anal(6)*h_mesh%gauss%rnorms(1,ls,ms))
4205  src_h(i,5) = src_h(i,5)+x*(jsolh_anal(3)*h_mesh%gauss%rnorms(1,ls,ms))
4206  src_h(i,6) = src_h(i,6)+x*(jsolh_anal(4)*h_mesh%gauss%rnorms(1,ls,ms))
4207 
4208  END DO
4209  ENDDO
4210 
4211  END DO
4212 
4213  IF (h_mesh%np /= 0) THEN
4214 !!$ ALLOCATE(idxn(H_mesh%np))
4215  idxn = la_h%loc_to_glob(1,:)-1
4216  CALL vecsetvalues(vb_1, h_mesh%np, idxn, src_h(:,1), add_values, ierr)
4217  CALL vecsetvalues(vb_2, h_mesh%np, idxn, src_h(:,2), add_values, ierr)
4218  idxn = la_h%loc_to_glob(2,:)-1
4219  CALL vecsetvalues(vb_1, h_mesh%np, idxn, src_h(:,4), add_values, ierr)
4220  CALL vecsetvalues(vb_2, h_mesh%np, idxn, src_h(:,3), add_values, ierr)
4221  idxn = la_h%loc_to_glob(3,:)-1
4222  CALL vecsetvalues(vb_1, h_mesh%np, idxn, src_h(:,5), add_values, ierr)
4223  CALL vecsetvalues(vb_2, h_mesh%np, idxn, src_h(:,6), add_values, ierr)
4224 !!$ DEALLOCATE(idxn)
4225  END IF
4226  CALL vecassemblybegin(vb_1,ierr)
4227  CALL vecassemblyend(vb_1,ierr)
4228  CALL vecassemblybegin(vb_2,ierr)
4229  CALL vecassemblyend(vb_2,ierr)
4230 
4231  END SUBROUTINE rhs_dirichlet
4232 
4233  SUBROUTINE dirichlet_cavities(communicator, interface_H_phi, mesh, js_D)
4234  USE def_type_mesh
4235  USE chaine_caractere
4236  USE input_data
4237  USE my_util
4238  IMPLICIT NONE
4239  TYPE(interface_type), INTENT(IN) :: interface_h_phi
4240  TYPE(mesh_type), INTENT(IN) :: mesh
4241  INTEGER, POINTER, DIMENSION(:) :: js_d
4242  INTEGER, ALLOCATABLE, DIMENSION(:) :: on_proc_loc, on_proc, not_cav_loc, not_cav
4243  INTEGER, ALLOCATABLE, DIMENSION(:) :: is_ok, j_tmp
4244  INTEGER, DIMENSION(1) :: loc
4245  INTEGER :: m, ms, i, nb_dom, idx, nb_cav, ni
4246  LOGICAL :: okay
4247 #include "petsc/finclude/petsc.h"
4248  mpi_comm, INTENT(IN) :: communicator
4249  petscint :: rank
4250  petscerrorcode :: ierr
4251 
4252  IF (inputs%nb_dom_phi==0) RETURN
4253 
4254  CALL mpi_comm_rank(communicator, rank, ierr)
4255 
4256  nb_dom = inputs%nb_dom_phi
4257  ALLOCATE(on_proc_loc(nb_dom), on_proc(nb_dom))
4258  ALLOCATE(not_cav_loc(nb_dom), not_cav(nb_dom))
4259  on_proc_loc = -1
4260  on_proc = -1
4261  not_cav_loc = -1
4262  not_cav = -1
4263 
4264  DO m = 1, mesh%me
4265  i = mesh%i_d(m)
4266  IF (minval(abs(inputs%list_dom_phi-i)) /= 0) THEN
4267  WRITE(*,*) 'error in dirichlet cavities'
4268  END IF
4269  loc = minloc(abs(inputs%list_dom_phi-i))
4270  on_proc_loc(loc(1)) = rank
4271  END DO
4272  IF (mesh%mes /= 0) THEN
4273  ALLOCATE(is_ok(mesh%mes))
4274  is_ok = mesh%i_d(mesh%neighs)
4275  IF (interface_h_phi%mes /=0) THEN
4276  is_ok(interface_h_phi%mesh2) = 0
4277  END IF
4278  DO ms = 1, mesh%mes
4279  IF (sum(abs(mesh%rr(1,mesh%jjs(:,ms)))) .LT. 1.d-12) THEN
4280  is_ok(ms) = 0
4281  END IF
4282  IF (inputs%my_periodic%nb_periodic_pairs /=0) THEN
4283  IF (minval(abs(inputs%my_periodic%list_periodic-mesh%sides(ms))) == 0) THEN
4284  is_ok(ms) = 0
4285  END IF
4286  END IF
4287  END DO
4288 
4289  DO ms = 1, mesh%mes
4290  IF (is_ok(ms) == 0) cycle
4291  i = is_ok(ms)
4292  IF (minval(abs(inputs%list_dom_phi-i)) /= 0) THEN
4293  WRITE(*,*) 'error in dirichlet cavities'
4294  END IF
4295  loc = minloc(abs(inputs%list_dom_phi-i))
4296  not_cav_loc(loc(1)) = rank
4297  END DO
4298  END IF
4299  CALL mpi_allreduce(on_proc_loc, on_proc, nb_dom, mpi_integer, mpi_max, communicator, ierr)
4300  CALL mpi_allreduce(not_cav_loc, not_cav, nb_dom, mpi_integer, mpi_max, communicator, ierr)
4301 
4302  ALLOCATE(j_tmp(SIZE(js_d)+nb_dom))
4303  j_tmp(1:SIZE(js_d)) = js_d
4304  idx = SIZE(js_d)
4305  DO i = 1, nb_dom
4306  IF ( (not_cav(i)==-1) .AND. (on_proc(i)==rank) ) THEN
4307  idx = idx + 1
4308  okay = .false.
4309  DO m = 1, mesh%me
4310  IF (mesh%i_d(m) == inputs%list_dom_phi(i)) THEN
4311  DO ni = 1, mesh%gauss%n_w
4312  IF (minval(abs(mesh%jjs-mesh%jj(ni,m))) /=0) THEN
4313  j_tmp(idx) = mesh%jj(ni,m)
4314  okay = .true.
4315  EXIT
4316  END IF
4317  END DO
4318  IF (okay) THEN
4319  WRITE(*,*) 'add ', j_tmp(idx), 'in dom ', inputs%list_dom_phi(i), ' : proc ', rank
4320  WRITE(*,*) 'add ', mesh%rr(:,j_tmp(idx)), mesh%i_d(m)
4321  EXIT
4322  END IF
4323  END IF
4324  END DO
4325  END IF
4326  END DO
4327 
4328  nb_cav = idx - SIZE(js_d)
4329  IF (nb_cav /= 0) THEN
4330  DEALLOCATE(js_d)
4331  ALLOCATE(js_d(idx))
4332  js_d = j_tmp(1:idx)
4333  END IF
4334 
4335  DEALLOCATE(on_proc_loc, on_proc, j_tmp)
4336  DEALLOCATE(not_cav_loc, not_cav)
4337  IF (ALLOCATED(is_ok)) DEALLOCATE(is_ok)
4338 
4339  WRITE(*,'(a,x,i2,x,a,x,i2)') 'I have detected', nb_cav, ' cavity(ies) on proc', rank
4340 
4341  END SUBROUTINE dirichlet_cavities
4342 
4343  SUBROUTINE smb_sigma_prod_curl(communicator, mesh, list_mode, H_in, sigma_bar, sigma_in, V_out)
4344  !=================================
4345  USE sft_parallele
4346  USE chaine_caractere
4347  USE input_data
4348  USE def_type_mesh
4349  USE user_data
4350  IMPLICIT NONE
4351 
4352  TYPE(mesh_type), INTENT(IN) :: mesh
4353  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
4354  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: h_in
4355  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: sigma_bar
4356  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: sigma_in
4357  REAL(KIND=8), DIMENSION(:,:,:) :: v_out
4358 
4359  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%me,6,SIZE(list_mode)) :: h_gauss, roth
4360  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%me,2,SIZE(list_mode)) :: sigma_gauss
4361  REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%me,6,SIZE(list_mode)) :: roth_bar
4362  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
4363  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
4364  INTEGER :: m, l , i, mode, index, k
4365  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,6) :: h_in_loc
4366  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,2) :: sigma_in_loc
4367  REAL(KIND=8) :: ray
4368  !REAL(KIND=8), DIMENSION(3) :: temps
4369  INTEGER :: nb_procs, bloc_size, m_max_pad, code
4370 #include "petsc/finclude/petsc.h"
4371  mpi_comm :: communicator
4372 
4373  DO i = 1, SIZE(list_mode)
4374  mode = list_mode(i)
4375  index = 0
4376  DO m = 1, mesh%me
4377  j_loc = mesh%jj(:,m)
4378  DO k = 1, 6
4379  h_in_loc(:,k) = h_in(j_loc,k,i)
4380  END DO
4381  DO k = 1, 2
4382  sigma_in_loc(:,k) = sigma_in(j_loc,k,i)
4383  END DO
4384 
4385  DO l = 1, mesh%gauss%l_G
4386  index = index + 1
4387  dw_loc = mesh%gauss%dw(:,:,l,m)
4388 
4389  !===Compute radius of Gauss point
4390  ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
4391 
4392  !-----------------magnetic field on gauss points---------------------------
4393  h_gauss(index,1,i) = sum(h_in_loc(:,1)*mesh%gauss%ww(:,l))
4394  h_gauss(index,3,i) = sum(h_in_loc(:,3)*mesh%gauss%ww(:,l))
4395  h_gauss(index,5,i) = sum(h_in_loc(:,5)*mesh%gauss%ww(:,l))
4396 
4397  h_gauss(index,2,i) = sum(h_in_loc(:,2)*mesh%gauss%ww(:,l))
4398  h_gauss(index,4,i) = sum(h_in_loc(:,4)*mesh%gauss%ww(:,l))
4399  h_gauss(index,6,i) = sum(h_in_loc(:,6)*mesh%gauss%ww(:,l))
4400  !-----------------sigma on gauss points------------------------------------
4401  sigma_gauss(index,1,i) = sum(sigma_in_loc(:,1)*mesh%gauss%ww(:,l))
4402  sigma_gauss(index,2,i) = sum(sigma_in_loc(:,2)*mesh%gauss%ww(:,l))
4403  !-----------------Curl of H on gauss points--------------------------------
4404  !coeff sur les cosinus
4405  roth(index,1,i) = mode/ray*h_gauss(index,6,i) &
4406  -sum(h_in_loc(:,3)*dw_loc(2,:))
4407  roth(index,4,i) = sum(h_in_loc(:,2)*dw_loc(2,:)) &
4408  -sum(h_in_loc(:,6)*dw_loc(1,:))
4409  roth(index,5,i) = 1/ray*h_gauss(index,3,i) &
4410  +sum(h_in_loc(:,3)*dw_loc(1,:)) &
4411  -mode/ray*h_gauss(index,2,i)
4412  !coeff sur les sinus
4413  roth(index,2,i) =-mode/ray*h_gauss(index,5,i) &
4414  -sum(h_in_loc(:,4)*dw_loc(2,:))
4415  roth(index,3,i) = sum(h_in_loc(:,1)*dw_loc(2,:)) &
4416  -sum(h_in_loc(:,5)*dw_loc(1,:))
4417  roth(index,6,i) = 1/ray*h_gauss(index,4,i) &
4418  +sum(h_in_loc(:,4)*dw_loc(1,:))&
4419  +mode/ray*h_gauss(index,1,i)
4420 
4421  DO k = 1, 6
4422  roth_bar(index,k,i) = roth(index,k,i)/sum(sigma_bar(j_loc)*mesh%gauss%ww(:,l))
4423  END DO
4424  ENDDO
4425  ENDDO
4426  END DO
4427 
4428  !temps = 0
4429  CALL mpi_comm_size(communicator, nb_procs, code)
4430  bloc_size = SIZE(sigma_gauss,1)/nb_procs+1
4431  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
4432 
4433  CALL fft_scalar_vect_dcl(communicator, roth, sigma_gauss, v_out, 2, nb_procs, &
4434  bloc_size, m_max_pad)
4435 
4436  v_out = roth_bar - v_out
4437 
4438  !write(*,*) ' Temps de Comm ', temps(1)
4439  !write(*,*) ' Temps de Calc ', temps(2)
4440  !write(*,*) ' Temps de Chan ', temps(3)
4441 
4442  END SUBROUTINE smb_sigma_prod_curl
4443 
4444  SUBROUTINE smb_sigma_prod_curl_bdy(communicator, mesh, Dirichlet_bdy_H_sides, list_mode, H_in, sigma_bar, sigma_in, V_out)
4445  !=================================
4446  USE sft_parallele
4447  USE chaine_caractere
4448  USE input_data
4449  USE def_type_mesh
4450  USE user_data
4451  IMPLICIT NONE
4452 
4453  TYPE(mesh_type), INTENT(IN) :: mesh
4454  INTEGER, DIMENSION(:), INTENT(IN) :: dirichlet_bdy_h_sides
4455  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
4456  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: h_in
4457  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: sigma_bar
4458  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: sigma_in
4459  REAL(KIND=8), DIMENSION(:,:,:) :: v_out
4460 
4461  REAL(KIND=8), DIMENSION(mesh%gauss%l_Gs*SIZE(Dirichlet_bdy_H_sides),6,SIZE(list_mode)) :: h_gauss, roth
4462  REAL(KIND=8), DIMENSION(mesh%gauss%l_Gs*SIZE(Dirichlet_bdy_H_sides),6,SIZE(list_mode)) :: roth_bar
4463  REAL(KIND=8), DIMENSION(mesh%gauss%l_Gs*SIZE(Dirichlet_bdy_H_sides),2,SIZE(list_mode)) :: sigma_gauss
4464  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
4465  INTEGER :: ms, ls , i, mode, index, k
4466  INTEGER, DIMENSION(mesh%gauss%n_ws) :: j_loc
4467  REAL(KIND=8), DIMENSION(mesh%gauss%n_ws,6) :: h_in_loc
4468  REAL(KIND=8), DIMENSION(mesh%gauss%n_ws,2) :: sigma_in_loc
4469  REAL(KIND=8) :: ray
4470  !REAL(KIND=8), DIMENSION(3) :: temps
4471  INTEGER :: nb_procs, bloc_size, m_max_pad, code, count, m1
4472 #include "petsc/finclude/petsc.h"
4473  mpi_comm :: communicator
4474 
4475  DO i = 1, SIZE(list_mode)
4476  mode = list_mode(i)
4477  index = 0
4478  DO count = 1, SIZE(dirichlet_bdy_h_sides)
4479  ms = dirichlet_bdy_h_sides(count)
4480  m1 = mesh%neighs(ms)
4481 
4482  j_loc = mesh%jjs(:,ms)
4483  DO k = 1, 6
4484  h_in_loc(:,k) = h_in(j_loc,k,i)
4485  END DO
4486  DO k = 1, 2
4487  sigma_in_loc(:,k) = sigma_in(j_loc,k,i)
4488  END DO
4489 
4490  DO ls = 1, mesh%gauss%l_Gs
4491  index = index + 1
4492  dw_loc = mesh%gauss%dw_s(:,:,ls,ms)
4493 
4494  !===Compute radius of Gauss point
4495  ray = sum(mesh%rr(1,mesh%jjs(:,ms))*mesh%gauss%wws(:,ls))
4496 
4497  !-----------------magnetic field on bdy gauss points---------------------------
4498  h_gauss(index,1,i) = sum(h_in_loc(:,1)*mesh%gauss%wws(:,ls))
4499  h_gauss(index,3,i) = sum(h_in_loc(:,3)*mesh%gauss%wws(:,ls))
4500  h_gauss(index,5,i) = sum(h_in_loc(:,5)*mesh%gauss%wws(:,ls))
4501 
4502  h_gauss(index,2,i) = sum(h_in_loc(:,2)*mesh%gauss%wws(:,ls))
4503  h_gauss(index,4,i) = sum(h_in_loc(:,4)*mesh%gauss%wws(:,ls))
4504  h_gauss(index,6,i) = sum(h_in_loc(:,6)*mesh%gauss%wws(:,ls))
4505  !-----------------sigma on bdy gauss points------------------------------------
4506  sigma_gauss(index,1,i) = sum(sigma_in_loc(:,1)*mesh%gauss%wws(:,ls))
4507  sigma_gauss(index,2,i) = sum(sigma_in_loc(:,2)*mesh%gauss%wws(:,ls))
4508  !-----------------Curl of H on bdy gauss points--------------------------------
4509  !coeff sur les cosinus
4510  roth(index,1,i) = mode/ray*h_gauss(index,6,i) &
4511 !!$ -SUM(H_in_loc(:,3)*dw_loc(2,:))
4512  -sum(h_in(mesh%jj(:,m1),3,i)*dw_loc(2,:))
4513 
4514 !!$ RotH(index,4,i) = SUM(H_in_loc(:,2)*dw_loc(2,:)) &
4515 !!$ -SUM(H_in_loc(:,6)*dw_loc(1,:))
4516  roth(index,4,i) = sum(h_in(mesh%jj(:,m1),2,i)*dw_loc(2,:)) &
4517  -sum(h_in(mesh%jj(:,m1),6,i)*dw_loc(1,:))
4518 
4519  roth(index,5,i) = 1/ray*h_gauss(index,3,i) &
4520 !!$ +SUM(H_in_loc(:,3)*dw_loc(1,:)) &
4521  +sum(h_in(mesh%jj(:,m1),3,i)*dw_loc(1,:)) &
4522  -mode/ray*h_gauss(index,2,i)
4523 
4524  !coeff sur les sinus
4525  roth(index,2,i) =-mode/ray*h_gauss(index,5,i) &
4526 !!$ -SUM(H_in_loc(:,4)*dw_loc(2,:))
4527  -sum(h_in(mesh%jj(:,m1),4,i)*dw_loc(2,:))
4528 
4529 !!$ RotH(index,3,i) = SUM(H_in_loc(:,1)*dw_loc(2,:)) &
4530 !!$ -SUM(H_in_loc(:,5)*dw_loc(1,:))
4531  roth(index,3,i) = sum(h_in(mesh%jj(:,m1),1,i)*dw_loc(2,:)) &
4532  -sum(h_in(mesh%jj(:,m1),5,i)*dw_loc(1,:))
4533 
4534  roth(index,6,i) = 1/ray*h_gauss(index,4,i) &
4535 !!$ +SUM(H_in_loc(:,4)*dw_loc(1,:))&
4536  +sum(h_in(mesh%jj(:,m1),4,i)*dw_loc(1,:))&
4537  +mode/ray*h_gauss(index,1,i)
4538 
4539  DO k = 1, 6
4540  roth_bar(index,k,i) = roth(index,k,i)/(sum(sigma_bar(j_loc)*mesh%gauss%wws(:,ls)))
4541  END DO
4542  END DO
4543  END DO
4544  END DO
4545 
4546  IF ( SIZE(dirichlet_bdy_h_sides).GE.1) THEN
4547  !temps = 0
4548  CALL mpi_comm_size(communicator, nb_procs, code)
4549  bloc_size = SIZE(sigma_gauss,1)/nb_procs+1
4550  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
4551 
4552  CALL fft_scalar_vect_dcl(communicator, roth, sigma_gauss, v_out, 2, nb_procs, &
4553  bloc_size, m_max_pad)
4554 
4555  v_out = roth_bar - v_out
4556  END IF
4557  !write(*,*) ' Temps de Comm ', temps(1)
4558  !write(*,*) ' Temps de Calc ', temps(2)
4559  !write(*,*) ' Temps de Chan ', temps(3)
4560 
4561  END SUBROUTINE smb_sigma_prod_curl_bdy
4562 
4563  SUBROUTINE smb_current_over_sigma(communicator, mesh, list_mode, &
4564  mu_h_field, mu_phi, sigma_tot, time, j_over_sigma)
4565  USE sft_parallele
4566  USE chaine_caractere
4567  USE input_data
4568  USE def_type_mesh
4569  USE boundary
4570  IMPLICIT NONE
4571  TYPE(mesh_type), INTENT(IN) :: mesh
4572  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
4573  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: sigma_tot
4574  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: mu_h_field
4575  REAL(KIND=8), INTENT(IN) :: mu_phi, time
4576  REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT) :: j_over_sigma
4577  REAL(KIND=8), DIMENSION(SIZE(mesh%rr,2),6,SIZE(list_mode)) :: j_node
4578  REAL(KIND=8), DIMENSION(6) :: b_ext_l
4579  REAL(KIND=8) :: muhl
4580  INTEGER :: mode, mesh_id1, k, i, n
4581  INTEGER :: nb_procs, bloc_size, m_max_pad, code
4582 #include "petsc/finclude/petsc.h"
4583  mpi_comm :: communicator
4584 
4585  !dummy values
4586  !ATTENTION the magnetic permeability needs to be constant for multiphase
4587  !and variable electrical condutivity
4588  !dummy values
4589  muhl=minval(mu_h_field)
4590  mesh_id1=1
4591  b_ext_l=1.d0
4592  !end dummy values
4593  DO i = 1, SIZE(list_mode)
4594  mode = list_mode(i)
4595  DO k = 1, 6
4596  DO n = 1, SIZE(mesh%rr,2)
4597  j_node(n,k,i) = jexact_gauss(k, mesh%rr(:,n), mode, mu_phi, 1.d0, muhl, &
4598  time, mesh_id1, b_ext_l)
4599  END DO
4600  END DO
4601  END DO
4602 
4603  CALL mpi_comm_size(communicator, nb_procs, code)
4604  bloc_size = SIZE(j_node,1)/nb_procs+1
4605  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
4606  CALL fft_scalar_vect_dcl(communicator, j_node, sigma_tot, j_over_sigma, 2, &
4607  nb_procs, bloc_size, m_max_pad)
4608 
4609  END SUBROUTINE smb_current_over_sigma
4610 
4611 END MODULE update_maxwell_with_h
4612 
def_type_mesh::interface_type
Definition: def_type_mesh.f90:113
user_data
Definition: read_user_data.f90:10
tn_axi
Definition: tn_axi.f90:5
update_maxwell_with_b::surf_int
subroutine surf_int(H_mesh, phi_mesh, interface_H_phi, interface_H_mu, list_dirichlet_sides_H, sigma, mu_phi, mu_H_field, time, mode, LA_H, LA_phi, vb_1, vb_2, R_fourier, index_fourier)
Definition: maxwell_update_time_with_B.f90:3601
solve_petsc
Definition: solver.f90:1
update_maxwell_with_b::courant_int_by_parts
subroutine courant_int_by_parts(H_mesh, phi_mesh, interface_H_phi, sigma, mu_phi, mu_H_field, time, mode, rhs_H, nl, LA_H, LA_phi, vb_1, vb_2, B_ext, H_pert, sigma_curl, J_over_sigma)
Definition: maxwell_update_time_with_B.f90:3146
sft_parallele
Definition: fft_parallel.f90:4
def_type_mesh::petsc_csr_la
Definition: def_type_mesh.f90:18
st_matrix::create_my_ghost
subroutine, public create_my_ghost(mesh, LA, ifrom)
Definition: st_csr.f90:14
dir_nodes
Definition: dir_nodes.f90:4
gauss_points
Definition: associate_gauss.f90:4
sub_plot
Definition: sub_plot.f90:4
dir_nodes_petsc::dirichlet_rhs
subroutine dirichlet_rhs(js_D, bs_D, b)
Definition: dir_nodes_petsc.f90:170
def_type_mesh::mesh_type
Definition: def_type_mesh.f90:73
st_matrix
Definition: st_csr.f90:4
boundary_test_21::sigma_bar_in_fourier_space
real(kind=8) function, dimension(size(h_mesh%rr, 2)), public sigma_bar_in_fourier_space(H_mesh)
Definition: condlim_test_21.f90:536
solve_petsc::init_solver
subroutine init_solver(my_par, my_ksp, matrix, communicator, solver, precond, opt_re_init)
Definition: solver.f90:11
update_maxwell_with_b::smb_sigma_prod_curl_bdy
subroutine smb_sigma_prod_curl_bdy(communicator, mesh, Dirichlet_bdy_H_sides, list_mode, H_in, sigma_bar, sigma_in, V_out)
Definition: maxwell_update_time_with_B.f90:5110
update_maxwell_with_b::mat_dirichlet_maxwell
subroutine mat_dirichlet_maxwell(H_mesh, Dirichlet_bdy_H_sides, mode, stab, mu_H_field, LA_H, H_p_phi_mat1, H_p_phi_mat2, sigma_np)
Definition: maxwell_update_time_with_B.f90:2692
dir_nodes_petsc::dirichlet_nodes_parallel
subroutine dirichlet_nodes_parallel(mesh, list_dirichlet_sides, js_d)
Definition: dir_nodes_petsc.f90:52
update_maxwell_with_b::mat_maxwell_mu
subroutine mat_maxwell_mu(H_mesh, interface_H_mu, mode, stab, mu_H_field, sigma, LA_H, H_p_phi_mat1, H_p_phi_mat2)
Definition: maxwell_update_time_with_B.f90:3886
dir_nodes_petsc
Definition: dir_nodes_petsc.f90:4
boundary_test_10::jexact_gauss
real(kind=8) function, public jexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, t, mesh_id, opt_B_ext)
Definition: condlim_test_10.f90:380
update_maxwell_with_h
Definition: maxwell_update_time_with_H.f90:6
solve_petsc::solver
subroutine solver(my_ksp, b, x, reinit, verbose)
Definition: solver.f90:95
def_type_mesh
Definition: def_type_mesh.f90:4
my_util
Definition: my_util.f90:1
dir_nodes_petsc::dirichlet_m_parallel
subroutine dirichlet_m_parallel(matrix, glob_js_D)
Definition: dir_nodes_petsc.f90:149
update_maxwell_with_h::maxwell_decouple_with_h
subroutine, public maxwell_decouple_with_h(comm_one_d, H_mesh, pmag_mesh, phi_mesh, interface_H_phi, interface_H_mu, Hn, Bn, phin, Hn1, Bn1, phin1, vel, stab_in, sigma_in, R_fourier, index_fourier, mu_H_field, mu_phi, time, dt, Rem, list_mode, H_phi_per, LA_H, LA_pmag, LA_phi, LA_mhd, sigma_ns_in, jj_v_to_H)
Definition: maxwell_update_time_with_H.f90:16
input_data
Definition: read_my_data.f90:200
def_type_mesh::periodic_type
Definition: def_type_mesh.f90:54
gauss_points::gauss
subroutine gauss(mesh)
Definition: associate_gauss.f90:23
chaine_caractere
Definition: chaine_caractere_par.f90:4
boundary_test_10::phiexact
real(kind=8) function, dimension(size(rr, 2)), public phiexact(TYPE, rr, m, mu_phi, t)
Definition: condlim_test_10.f90:360
update_maxwell_with_b::mat_h_p_phi_maxwell
subroutine mat_h_p_phi_maxwell(H_mesh, pmag_mesh, phi_mesh, interface_H_phi, mode, mu_H_field, mu_phi, c_mass, stab, R_fourier, index_fourier, LA_H, LA_pmag, LA_phi, H_p_phi_mat1, H_p_phi_mat2, sigma_np)
Definition: maxwell_update_time_with_B.f90:713
sft_parallele_obsolete::ref
subroutine, public ref(communicator, V1_in, V2_in, V_out, temps)
Definition: fft_parallel_obsolete.f90:2259
sft_parallele::fft_scalar_vect_dcl
subroutine, public fft_scalar_vect_dcl(communicator, V1_in, V2_in, V_out, pb, nb_procs, bloc_size, m_max_pad, temps)
Definition: fft_parallel.f90:853
tn_axi::norm_sf
real(kind=8) function norm_sf(communicator, norm_type, mesh, list_mode, v)
Definition: tn_axi.f90:38
dir_nodes::dirichlet_nodes
subroutine dirichlet_nodes(jjs_in, sides_in, dir_in, js_d)
Definition: dir_nodes.f90:495
prep_maill
Definition: prep_mesh.f90:5
boundary_test_10::eexact_gauss
real(kind=8) function, public eexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, t)
Definition: condlim_test_10.f90:402
my_util::user_time
real(kind=8) function user_time()
Definition: my_util.f90:7
periodic::periodic_rhs_petsc
subroutine, public periodic_rhs_petsc(n_bord, list, perlist, v_rhs, LA)
Definition: prep_mesh_periodic.f90:455
update_maxwell_with_b::courant_mu
subroutine courant_mu(H_mesh, interface_H_mu, sigma, mu_H_field, time, mode, nl, LA_H, vb_1, vb_2, B_ext)
Definition: maxwell_update_time_with_B.f90:4553
update_maxwell_with_b::rhs_dirichlet
subroutine rhs_dirichlet(H_mesh, Dirichlet_bdy_H_sides, sigma, mu_H_field, time, mode, nl, stab, LA_H, vb_1, vb_2, B_ext, J_over_sigma, sigma_curl_bdy_in)
Definition: maxwell_update_time_with_B.f90:4722
update_maxwell_with_b::dirichlet_cavities
subroutine dirichlet_cavities(communicator, interface_H_phi, mesh, js_D)
Definition: maxwell_update_time_with_B.f90:4875
my_util::error_petsc
subroutine error_petsc(string)
Definition: my_util.f90:15
sft_parallele::fft_par_cross_prod_dcl
subroutine, public fft_par_cross_prod_dcl(communicator, V1_in, V2_in, V_out, nb_procs, bloc_size, m_max_pad, temps)
Definition: fft_parallel.f90:144
periodic::periodic_matrix_petsc
subroutine, public periodic_matrix_petsc(n_bord, list, perlist, matrix, LA)
Definition: prep_mesh_periodic.f90:377
st_matrix::extract
subroutine, public extract(xghost, ks, ke, LA, phi)
Definition: st_csr.f90:33
solve_petsc::create_local_petsc_matrix
subroutine create_local_petsc_matrix(communicator, LA, matrix, clean)
Definition: solver.f90:142
boundary_test_10::hexact
real(kind=8) function, dimension(size(rr, 2)), public hexact(H_mesh, TYPE, rr, m, mu_H_field, t)
Definition: condlim_test_10.f90:316
z
section doc_intro_frame_work_num_app Numerical approximation subsection doc_intro_fram_work_num_app_Fourier_FEM Fourier Finite element representation The SFEMaNS code uses a hybrid Fourier Finite element formulation The Fourier decomposition allows to approximate the problem’s solutions for each Fourier mode modulo nonlinear terms that are made explicit The variables are then approximated on a meridian section of the domain with a finite element method The numerical approximation of a function f $f f is written in the following generic z
Definition: doc_intro.h:193
update_maxwell_with_b::smb_current_over_sigma
subroutine smb_current_over_sigma(communicator, mesh, list_mode, mu_H_field, mu_phi, sigma_tot, time, J_over_sigma)
Definition: maxwell_update_time_with_B.f90:5222
update_maxwell_with_b::smb_sigma_prod_curl
subroutine smb_sigma_prod_curl(communicator, mesh, list_mode, H_in, sigma_bar, sigma_in, V_out)
Definition: maxwell_update_time_with_B.f90:5017