SFEMaNS  version 5.3
Reference documentation for SFEMaNS
metis_distribution_sfemans.f90
Go to the documentation of this file.
1 MODULE metis_sfemans
2 #include "petsc/finclude/petsc.h"
3  USE petsc
4  PUBLIC :: part_mesh_mhd, free_mesh, extract_mesh, free_interface, part_mesh_mhd_bis, part_mesh_m_t_h_phi
5  PRIVATE
6  REAL(KIND=8) :: epsilon = 1.d-10
7 !!$ Dummy for metis...
8  INTEGER :: metis_noptions=40, metis_option_numbering=18
9  !==Go see in metis.h the actual values assigned to METIS_NOPTIONS, METIS_OPTION_NUMBERING
10 !!$ Dummy for metis...
11 CONTAINS
12  !===Convention: Domain NS \subset domain Temp \subset Domain H
13  SUBROUTINE part_mesh_m_t_h_phi(nb_proc,list_u, list_T_in, list_h_in, list_phi ,mesh,list_of_interfaces,part,my_periodic)
14  USE def_type_mesh
15  USE my_util
16  USE sub_plot
17  USE periodic
18  IMPLICIT NONE
19  TYPE(mesh_type) :: mesh
20  INTEGER, DIMENSION(mesh%me) :: part
21  INTEGER, DIMENSION(:) :: list_of_interfaces
22  INTEGER, DIMENSION(:) :: list_u, list_T_in, list_h_in, list_phi
23  TYPE(periodic_data), OPTIONAL :: my_periodic
24 
25  LOGICAL, DIMENSION(mesh%mes) :: virgins
26  INTEGER, DIMENSION(3,mesh%me) :: neigh_new
27  INTEGER, DIMENSION(5) :: opts
28  INTEGER, DIMENSION(SIZE(mesh%jjs,1)) :: i_loc
29  INTEGER, DIMENSION(:), ALLOCATABLE :: xadj_u, xadj_T, xadj_h, xadj_phi, list_h, list_T
30  INTEGER, DIMENSION(:), ALLOCATABLE :: xind_u, xind_T, xind_h, xind_phi
31  INTEGER, DIMENSION(:), ALLOCATABLE :: vwgt, adjwgt
32  INTEGER, DIMENSION(:), ALLOCATABLE :: u2glob, T2glob, h2glob, phi2glob
33  INTEGER, DIMENSION(:), ALLOCATABLE :: part_u, part_T, part_h, part_phi
34  INTEGER, DIMENSION(1) :: jm_loc
35  INTEGER, DIMENSION(mesh%np,3) :: per_pts
36  INTEGER, DIMENSION(mesh%me) :: glob2loc
37  INTEGER, DIMENSION(mesh%np) :: indicator
38  INTEGER, DIMENSION(3) :: j_loc
39  INTEGER :: nb_neigh, edge, m, ms, n, nb, numflag, p, wgtflag, j, &
40  ns, nws, msop, nsop, proc, iop, mop, s2, k, me_u, me_T, me_h, me_phi, idm
41  REAL(KIND=8) :: err
42  LOGICAL :: test
43  !===(JLG) Feb 20, 2019. Petsc developpers decide to use REAL(KIND=4) to interface with metis
44  !REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: tpwgts
45  !REAL(KIND=8), DIMENSION(1) :: ubvec
46  REAL(KIND=4), DIMENSION(:), ALLOCATABLE :: tpwgts
47  REAL(KIND=4), DIMENSION(1) :: ubvec
48  REAL(KIND=4) :: one_K4=1.0
49  !===(JLG)Feb 20, 2019.
50  INTEGER, DIMENSION(METIS_NOPTIONS) :: metis_opt
51  petscmpiint :: nb_proc
52 !!$ WARNING, FL 1/2/13 : TO BE ADDED IN NEEDED
53  petscerrorcode :: ierr
54  petscmpiint :: rank
55  CALL mpi_comm_rank(mpi_comm_world,rank,ierr)
56 !!$ WARNING, FL 1/2/13 : TO BE ADDED IN NEEDED
57 
58  IF (nb_proc==1) THEN
59  part = 1
60  RETURN
61  END IF
62  glob2loc = 0
63 
64  !===Create list_T = list_T_in \ list_u
65  nb = 0
66  DO j = 1, SIZE(list_t_in)
67  IF (minval(abs(list_t_in(j)-list_u))==0) cycle
68  nb = nb + 1
69  END DO
70  ALLOCATE(list_t(nb))
71  nb = 0
72  DO j = 1, SIZE(list_t_in)
73  IF (minval(abs(list_t_in(j)-list_u))==0) cycle
74  nb = nb + 1
75  list_t(nb) = list_t_in(j)
76  END DO
77  !==Done with list_T
78 
79  !===Create list_h = list_h_in \ list_T_in
80  nb = 0
81  DO j = 1, SIZE(list_h_in)
82  IF (minval(abs(list_h_in(j)-list_t_in))==0) cycle
83  IF (minval(abs(list_h_in(j)-list_u))==0) cycle
84  nb = nb + 1
85  END DO
86  ALLOCATE(list_h(nb))
87  nb = 0
88  DO j = 1, SIZE(list_h_in)
89  IF (minval(abs(list_h_in(j)-list_t_in))==0) cycle
90  IF (minval(abs(list_h_in(j)-list_u))==0) cycle
91  nb = nb + 1
92  list_h(nb) = list_h_in(j)
93  END DO
94  !==Done with list_h
95 
96  !===Create neigh_new for interfaces
97  nws = SIZE( mesh%jjs,1)
98  neigh_new = mesh%neigh
99  IF (SIZE(list_of_interfaces)/=0) THEN
100  virgins = .true.
101  DO ms = 1, mesh%mes
102  IF (.NOT.virgins(ms)) cycle
103  IF (minval(abs(mesh%sides(ms)-list_of_interfaces))/=0) cycle !==ms not on a cut
104  i_loc = mesh%jjs(:,ms)
105  DO msop = 1, mesh%mes
106  IF (msop==ms .OR. .NOT.virgins(msop)) cycle
107  IF (minval(abs(mesh%sides(msop)-list_of_interfaces))/=0) cycle !==msop not on a cut
108  DO ns = 1, nws
109  test = .false.
110  DO nsop = 1, nws
111  iop = mesh%jjs(nsop,msop)
112  IF (maxval(abs(mesh%rr(:,i_loc(ns))-mesh%rr(:,iop))).LT.epsilon) THEN
113  test = .true.
114  EXIT
115  END IF
116  END DO
117  IF (.NOT.test) THEN
118  EXIT !==This msop does not coincide with ms
119  END IF
120  END DO
121  IF (test) EXIT
122  END DO
123  IF (.NOT.test) THEN
124  CALL error_petsc(.NOT.'BUG in part_mesh_M_T_H_phi, test ')
125  END IF
126  DO n = 1, 3
127  IF (neigh_new(n,mesh%neighs(msop))==0) THEN
128  neigh_new(n,mesh%neighs(msop)) = mesh%neighs(ms)
129  END IF
130  IF (neigh_new(n,mesh%neighs(ms))==0) THEN
131  neigh_new(n,mesh%neighs(ms)) = mesh%neighs(msop)
132  END IF
133  END DO
134  virgins(ms) = .false.
135  virgins(msop) = .false.
136  END DO
137  END IF
138  !===End Create neigh_new for interfaces
139 
140  !===Create neigh_new for periodic faces
141  IF (PRESENT(my_periodic)) THEN
142  IF (my_periodic%nb_periodic_pairs/=0) THEN
143  DO ms = 1, mesh%mes
144  m = mesh%neighs(ms)
145  IF (minval(abs(mesh%sides(ms)-my_periodic%list_periodic(1,:))) == 0) THEN
146  jm_loc = minloc(abs(mesh%sides(ms)-my_periodic%list_periodic(1,:)))
147  s2 = my_periodic%list_periodic(2,jm_loc(1))
148  test = .false.
149  DO msop = 1, mesh%mes
150  IF (mesh%sides(msop) /= s2) cycle
151 
152  err = 0.d0
153  DO ns = 1, SIZE(my_periodic%vect_e,1)
154  err = err + abs(sum(mesh%rr(ns,mesh%jjs(:,ms))-mesh%rr(ns,mesh%jjs(:,msop)) &
155  +my_periodic%vect_e(ns,jm_loc(1))))
156  END DO
157 
158  IF (err .LE. epsilon) THEN
159  test = .true.
160  EXIT
161  END IF
162  END DO
163  IF (.NOT.test) THEN
164  CALL error_petsc('BUG in part_mesh_M_T_H_phi, mop not found')
165  END IF
166  mop = mesh%neighs(msop)
167  DO n = 1, 3
168  IF (neigh_new(n,m) == 0) THEN
169  neigh_new(n,m) = mop
170  END IF
171  IF (neigh_new(n,mop) == 0) THEN
172  neigh_new(n,mop) = m
173  END IF
174  END DO
175  END IF
176  END DO
177  END IF
178  END IF
179  !===End Create neigh_new for periodic faces
180 
181  !===Create glob2loc and u2glob, T2glob, h2glob, phi2glob
182  me_u = 0
183  me_t = 0
184  me_h = 0
185  me_phi = 0
186  DO m = 1, mesh%me
187  idm = mesh%i_d(m)
188  IF (minval(abs(idm-list_u))==0) THEN
189  me_u = me_u + 1
190  ELSE IF (minval(abs(idm-list_t))==0) THEN
191  me_t = me_t + 1
192  ELSE IF (minval(abs(idm-list_h))==0) THEN
193  me_h = me_h + 1
194  ELSE IF (minval(abs(idm-list_phi))==0) THEN
195  me_phi = me_phi + 1
196  ELSE
197  CALL error_petsc('BUG in part_mesh_M_T_H_phi : element not in the mesh')
198  END IF
199  END DO
200  ALLOCATE(u2glob(me_u), t2glob(me_t), h2glob(me_h), phi2glob(me_phi))
201  me_u = 0
202  me_t = 0
203  me_h = 0
204  me_phi = 0
205  DO m = 1, mesh%me
206  idm = mesh%i_d(m)
207  IF (minval(abs(idm-list_u))==0) THEN
208  me_u = me_u + 1
209  u2glob(me_u) = m
210  glob2loc(m) = me_u
211  ELSE IF (minval(abs(idm-list_t))==0) THEN
212  me_t = me_t + 1
213  t2glob(me_t) = m
214  glob2loc(m) = me_t
215  ELSE IF (minval(abs(idm-list_h))==0) THEN
216  me_h = me_h + 1
217  h2glob(me_h) = m
218  glob2loc(m) = me_h
219  ELSE IF (minval(abs(idm-list_phi))==0) THEN
220  me_phi = me_phi + 1
221  phi2glob(me_phi) = m
222  glob2loc(m) = me_phi
223  ELSE
224  CALL error_petsc('BUG in part_mesh_M_T_H_phi: element not in the mesh')
225  END IF
226  END DO
227  !===End Create glob2loc and u2glob, T2glob, h2glob, phi2glob
228 
229  !===Create the connectivity arrays Xind and Xadj based on neigh (for Metis)
230  nb_neigh = SIZE(mesh%neigh,1)
231  ALLOCATE(xind_u(me_u+1), xind_t(me_t+1), xind_h(me_h+1), xind_phi(me_phi+1))
232  xind_u(1) = 1
233  DO k = 1, me_u
234  m = u2glob(k)
235  nb = 0
236  DO n = 1, nb_neigh
237  mop = neigh_new(n,m)
238  IF (mop==0) cycle
239  IF (minval(abs(mesh%i_d(mop)-list_u))/=0) cycle
240  nb = nb + 1
241  END DO
242  xind_u(k+1) = xind_u(k) + nb
243  END DO
244  xind_t(1) = 1
245  DO k = 1, me_t
246  m = t2glob(k)
247  nb = 0
248  DO n = 1, nb_neigh
249  mop = neigh_new(n,m)
250  IF (mop==0) cycle
251  IF (minval(abs(mesh%i_d(mop)-list_t))/=0) cycle
252  nb = nb + 1
253  END DO
254  xind_t(k+1) = xind_t(k) + nb
255  END DO
256  xind_h(1) = 1
257  DO k = 1, me_h
258  m = h2glob(k)
259  nb = 0
260  DO n = 1, nb_neigh
261  mop = neigh_new(n,m)
262  IF (mop==0) cycle
263  IF (minval(abs(mesh%i_d(mop)-list_h))/=0) cycle
264  nb = nb + 1
265  END DO
266  xind_h(k+1) = xind_h(k) + nb
267  END DO
268  xind_phi(1) = 1
269  DO k = 1, me_phi
270  m = phi2glob(k)
271  nb = 0
272  DO n = 1, nb_neigh
273  mop = neigh_new(n,m)
274  IF (mop==0) cycle
275  IF (minval(abs(mesh%i_d(mop)-list_phi))/=0) cycle
276  nb = nb + 1
277  END DO
278  xind_phi(k+1) = xind_phi(k) + nb
279  END DO
280 
281  ALLOCATE(xadj_u(xind_u(me_u+1)-1))
282  ALLOCATE(xadj_t(xind_t(me_t+1)-1))
283  ALLOCATE(xadj_h(xind_h(me_h+1)-1))
284  ALLOCATE(xadj_phi(xind_phi(me_phi+1)-1))
285  p = 0
286  DO k = 1, me_u
287  m = u2glob(k)
288  DO n = 1, nb_neigh
289  mop = neigh_new(n,m)
290  IF (mop==0) cycle
291  IF (minval(abs(mesh%i_d(mop)-list_u))/=0) cycle
292  p = p + 1
293  xadj_u(p) = glob2loc(mop)
294  END DO
295  END DO
296  IF (p/=xind_u(me_u+1)-1) THEN
297  CALL error_petsc('BUG in part_mesh_M_T_H_phi, p/=xind_u(me_u+1)-1')
298  END IF
299  p = 0
300  DO k = 1, me_t
301  m = t2glob(k)
302  DO n = 1, nb_neigh
303  mop = neigh_new(n,m)
304  IF (mop==0) cycle
305  IF (minval(abs(mesh%i_d(mop)-list_t))/=0) cycle
306  p = p + 1
307  xadj_t(p) = glob2loc(mop)
308  END DO
309  END DO
310  IF (p/=xind_t(me_t+1)-1) THEN
311  CALL error_petsc('BUG in part_mesh_M_T_H_phi, p/=xind_T(me_T+1)-1')
312  END IF
313  p = 0
314  DO k = 1, me_h
315  m = h2glob(k)
316  DO n = 1, nb_neigh
317  mop = neigh_new(n,m)
318  IF (mop==0) cycle
319  IF (minval(abs(mesh%i_d(mop)-list_h))/=0) cycle
320  p = p + 1
321  xadj_h(p) = glob2loc(mop)
322  END DO
323  END DO
324  IF (p/=xind_h(me_h+1)-1) THEN
325  CALL error_petsc('BUG in part_mesh_M_T_H_phi, p/=xind_h(me_h+1)-1')
326  END IF
327  p = 0
328  DO k = 1, me_phi
329  m = phi2glob(k)
330  DO n = 1, nb_neigh
331  mop = neigh_new(n,m)
332  IF (mop==0) cycle
333  IF (minval(abs(mesh%i_d(mop)-list_phi))/=0) cycle
334  p = p + 1
335  xadj_phi(p) = glob2loc(mop)
336  END DO
337  END DO
338  IF (p/=xind_phi(me_phi+1)-1) THEN
339  CALL error_petsc('BUG in part_mesh_M_T_H_phi, p/=xind_phi(me_phi+1)-1')
340  END IF
341  !===End Create the connectivity arrays Xind and Xadj based on neigh (for Metis)
342 
343  !===Create partitions
344  opts = 0
345  numflag = 1
346  wgtflag = 2
347  ALLOCATE(tpwgts(nb_proc))
348  tpwgts=one_k4/nb_proc
349  CALL metis_setdefaultoptions(metis_opt)
350  metis_opt(metis_option_numbering)=1
351  ubvec=1.01
352  IF (me_u /= 0) THEN
353  ALLOCATE(vwgt(me_u), adjwgt(SIZE(xadj_u)), part_u(me_u))
354  vwgt = 1
355  adjwgt = 1
356  CALL metis_partgraphrecursive(me_u, 1, xind_u, xadj_u, vwgt, vwgt, adjwgt, nb_proc, tpwgts, ubvec, metis_opt, edge, part_u)
357  END IF
358  IF (me_t /= 0) THEN
359  IF (ALLOCATED(vwgt)) THEN
360  DEALLOCATE(vwgt, adjwgt)
361  END IF
362  ALLOCATE(vwgt(me_t), adjwgt(SIZE(xadj_t)), part_t(me_t))
363  vwgt = 1
364  adjwgt = 1
365  CALL metis_partgraphrecursive(me_t, 1, xind_t, xadj_t, vwgt, vwgt, adjwgt, nb_proc, tpwgts, ubvec, metis_opt, edge, part_t)
366  END IF
367  IF (me_h /= 0) THEN
368  IF (ALLOCATED(vwgt)) THEN
369  DEALLOCATE(vwgt, adjwgt)
370  END IF
371  ALLOCATE(vwgt(me_h), adjwgt(SIZE(xadj_h)), part_h(me_h))
372  vwgt = 1
373  adjwgt = 1
374  CALL metis_partgraphrecursive(me_h, 1, xind_h, xadj_h, vwgt, vwgt, adjwgt, nb_proc,tpwgts, ubvec, metis_opt, edge, part_h)
375  END IF
376  IF (me_phi /= 0) THEN
377  IF (ALLOCATED(vwgt)) THEN
378  DEALLOCATE(vwgt, adjwgt)
379  END IF
380  ALLOCATE(vwgt(me_phi), adjwgt(SIZE(xadj_phi)), part_phi(me_phi))
381  vwgt = 1
382  adjwgt = 1
383  CALL metis_partgraphrecursive(me_phi, 1, xind_phi,xadj_phi,vwgt, vwgt, adjwgt, nb_proc,tpwgts, &
384  ubvec, metis_opt, edge, part_phi)
385  END IF
386  !===Create global partition 'part'
387  part = -1
388  IF (me_u/=0) THEN
389  part(u2glob(:)) = part_u
390  END IF
391  IF (me_t/=0) THEN
392  part(t2glob(:)) = part_t
393  END IF
394  IF (me_h/=0) THEN
395  part(h2glob(:)) = part_h
396  END IF
397  IF (me_phi/=0) THEN
398  part(phi2glob(:)) = part_phi
399  END IF
400  IF (minval(part)==-1) THEN
401  CALL error_petsc('BUG in part_mesh_mhd_bis, MINVAL(part) == -1')
402  END IF
403  !===End Create global partition
404 
405  !===Create parts and modify part
406  !===Search on the boundary whether ms is on a cut.
407  IF (SIZE(mesh%jj,1)/=3) THEN
408  write(*,*) 'SIZE(mesh%jj,1)', SIZE(mesh%jj,1)
409  CALL error_petsc('BUG in part_mesh_M_T_H_phi, SIZE(mesh%jj,1)/=3')
410  END IF
411  indicator = -1
412  nws = SIZE( mesh%jjs,1)
413  IF (SIZE(list_of_interfaces)/=0) THEN
414  virgins = .true.
415  DO ms = 1, mesh%mes
416  IF (.NOT.virgins(ms)) cycle
417  IF (minval(abs(mesh%sides(ms)-list_of_interfaces))/=0) cycle !==ms not on a cut
418  i_loc = mesh%jjs(:,ms)
419  DO msop = 1, mesh%mes
420  IF (msop==ms .OR. .NOT.virgins(msop)) cycle
421  IF (minval(abs(mesh%sides(msop)-list_of_interfaces))/=0) cycle !==msop not on a cut
422  DO ns = 1, nws
423  test = .false.
424  DO nsop = 1, nws
425  iop = mesh%jjs(nsop,msop)
426  IF (maxval(abs(mesh%rr(:,i_loc(ns))-mesh%rr(:,iop))).LT.epsilon) THEN
427  test = .true.
428  EXIT
429  END IF
430  END DO
431  IF (.NOT.test) THEN
432  EXIT !==This msop does not coincide with ms
433  END IF
434  END DO
435  IF (test) EXIT
436  END DO
437  IF (.NOT.test) THEN
438  CALL error_petsc(.NOT.'BUG in part_mesh_M_T_H_phi, test ')
439  END IF
440  IF (part(mesh%neighs(ms)) == part(mesh%neighs(msop))) cycle !==ms is an internal cut
441  proc = min(part(mesh%neighs(ms)),part(mesh%neighs(msop)))
442  part(mesh%neighs(ms)) = proc !make sure interface are internal
443  part(mesh%neighs(msop)) = proc !make sure interface are internal
444  virgins(ms) = .false.
445  virgins(msop) = .false.
446  indicator(mesh%jjs(:,ms)) = proc
447  indicator(mesh%jjs(:,msop)) = proc
448  END DO
449  END IF
450  !===Fix the partition so that all the cells having one vertex on an
451  !===interface belong to the same processor as those sharing this vertices and
452  !===having two vertices on the interface (JLG + DCQ July 22 2015)
453  DO m = 1, mesh%me
454  j_loc = mesh%jj(:,m)
455  n = maxval(indicator(j_loc))
456  IF (n == -1) cycle
457  IF (indicator(j_loc(1))*indicator(j_loc(2))*indicator(j_loc(3))<0) cycle
458  part(m) = n
459  END DO
460  !===End create parts and modify part
461 
462  !===Move the two elements with one periodic face on same processor
463  IF (PRESENT(my_periodic)) THEN
464  IF (my_periodic%nb_periodic_pairs/=0) THEN
465  DO j = 1, mesh%np
466  per_pts(j,1) = j
467  END DO
468  per_pts(:,2:3) = 0
469  DO ms = 1, mesh%mes
470  m = mesh%neighs(ms)
471  IF ((minval(abs(mesh%sides(ms)-my_periodic%list_periodic(1,:))) /=0) .AND. &
472  (minval(abs(mesh%sides(ms)-my_periodic%list_periodic(2,:))) /=0) ) cycle
473  DO ns = 1, SIZE(mesh%jjs,1)
474  j = mesh%jjs(ns,ms)
475  per_pts(j,2) = m
476  DO msop = 1, mesh%mes
477  IF (minval(abs(mesh%sides(msop)-my_periodic%list_periodic(:,:))) /=0 ) cycle
478  IF (msop == ms) cycle
479  DO nsop = 1, SIZE(mesh%jjs,1)
480  IF (mesh%jjs(nsop,msop)==j) THEN
481  per_pts(j,3) = mesh%neighs(msop)
482  END IF
483  END DO
484  END DO
485  END DO
486  END DO
487  CALL reassign_per_pts(mesh, part, per_pts)
488  DO ms = 1, mesh%mes
489  m = mesh%neighs(ms)
490  IF (minval(abs(mesh%sides(ms)-my_periodic%list_periodic(1,:))) /= 0) cycle
491  jm_loc = minloc(abs(mesh%sides(ms)-my_periodic%list_periodic(1,:)))
492  s2 = my_periodic%list_periodic(2,jm_loc(1))
493  test = .false.
494  DO msop = 1, mesh%mes
495  IF (mesh%sides(msop) /= s2) cycle
496  err = 0.d0
497  DO ns = 1, SIZE(my_periodic%vect_e,1)
498  err = err + abs(sum(mesh%rr(ns,mesh%jjs(:,ms))-mesh%rr(ns,mesh%jjs(:,msop)) &
499  +my_periodic%vect_e(ns,jm_loc(1))))
500  END DO
501  IF (err .LE. epsilon) THEN
502  test = .true.
503  EXIT
504  END IF
505  END DO
506  IF (.NOT.test) THEN
507  CALL error_petsc('BUG in part_mesh_M_T_H_phi, mop not found')
508  END IF
509  IF (part(mesh%neighs(ms)) /= part(mesh%neighs(msop))) THEN !==ms is an internal cut
510  proc = min(part(mesh%neighs(ms)),part(mesh%neighs(msop)))
511  part(mesh%neighs(ms)) = proc !make sure interface are internal
512  part(mesh%neighs(msop)) = proc !make sure interface are internal
513  END IF
514  END DO
515  END IF
516  END IF
517  !===End Move the two elements with one periodic face on same processor
518 
519 !!$ WARNING, FL 1/2/13 : TO BE ADDED IF NEEDED
520  !================================================
521  IF (rank==0) THEN
522  CALL plot_const_p1_label(mesh%jj, mesh%rr, 1.d0*part, 'dd.plt')
523  END IF
524  !================================================
525 !!$ WARNING, FL 1/2/13 : TO BE ADDED IF NEEDED
526 
527  DEALLOCATE(vwgt,adjwgt)
528  IF (ALLOCATED(xadj_u)) DEALLOCATE(xadj_u)
529  IF (ALLOCATED(xadj_t)) DEALLOCATE(xadj_t)
530  IF (ALLOCATED(xadj_h)) DEALLOCATE(xadj_h)
531  IF (ALLOCATED(xadj_phi)) DEALLOCATE(xadj_phi)
532  IF (ALLOCATED(list_t)) DEALLOCATE(list_t)
533  IF (ALLOCATED(list_h)) DEALLOCATE(list_h)
534  IF (ALLOCATED(xind_u)) DEALLOCATE(xind_u)
535  IF (ALLOCATED(xind_t)) DEALLOCATE(xind_t)
536  IF (ALLOCATED(xind_h)) DEALLOCATE(xind_h)
537  IF (ALLOCATED(xind_phi)) DEALLOCATE(xind_phi)
538  IF (ALLOCATED(u2glob)) DEALLOCATE(u2glob)
539  IF (ALLOCATED(t2glob)) DEALLOCATE(t2glob)
540  IF (ALLOCATED(h2glob)) DEALLOCATE(h2glob)
541  IF (ALLOCATED(phi2glob)) DEALLOCATE(phi2glob)
542  IF (ALLOCATED(part_u)) DEALLOCATE(part_u)
543  IF (ALLOCATED(part_t)) DEALLOCATE(part_t)
544  IF (ALLOCATED(part_h)) DEALLOCATE(part_h)
545  IF (ALLOCATED(part_phi)) DEALLOCATE(part_phi)
546 
547  DEALLOCATE(tpwgts)
548 
549  END SUBROUTINE part_mesh_m_t_h_phi
550 
551  SUBROUTINE part_mesh_mhd(nb_proc,vwgt,mesh,list_of_interfaces,part,my_periodic)
552  USE def_type_mesh
553  USE my_util
554  USE sub_plot
555  USE periodic
556  IMPLICIT NONE
557  TYPE(mesh_type) :: mesh
558  INTEGER, DIMENSION(mesh%me+1) :: xind
559  INTEGER, DIMENSION(mesh%me) :: vwgt, part
560  INTEGER, DIMENSION(:) :: list_of_interfaces
561  TYPE(periodic_data), OPTIONAL :: my_periodic
562 
563  LOGICAL, DIMENSION(mesh%mes) :: virgins
564  INTEGER, DIMENSION(3,mesh%me) :: neigh_new
565  INTEGER, DIMENSION(5) :: opts
566  INTEGER, DIMENSION(SIZE(mesh%jjs,1)) :: i_loc
567  INTEGER, DIMENSION(:), ALLOCATABLE :: xadj, adjwgt
568  INTEGER, DIMENSION(1) :: jm_loc
569  INTEGER, DIMENSION(mesh%np,3) :: per_pts
570  INTEGER :: nb_neigh, edge, m, ms, n, nb, numflag, p, wgtflag, j, &
571  ns, nws, msop, nsop, proc, iop, mop, s2
572  REAL(KIND=8) :: err
573  !===(JLG) Feb 20, 2019. Petsc developpers decide to use REAL(KIND=4) to interface with metis
574  !REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: tpwgts
575  !REAL(KIND=8), DIMENSION(1) :: ubvec
576  REAL(KIND=4), DIMENSION(:), ALLOCATABLE :: tpwgts
577  REAL(KIND=4), DIMENSION(1) :: ubvec
578  REAL(KIND=4) :: one_K4=1.0
579  !===(JLG)Feb 20, 2019.
580  INTEGER, DIMENSION(METIS_NOPTIONS) :: metis_opt
581  LOGICAL :: test
582  petscmpiint :: nb_proc
583 !!$ WARNING, FL 1/2/13 : TO BE ADDED IN NEEDED
584  petscerrorcode :: ierr
585  petscmpiint :: rank
586 !!$ WARNING, FL 1/2/13 : TO BE ADDED IN NEEDED
587 
588  IF (nb_proc==1) THEN
589  part = 1
590  RETURN
591  END IF
592  nws = SIZE( mesh%jjs,1)
593  neigh_new = mesh%neigh
594  IF (SIZE(list_of_interfaces)/=0) THEN
595  virgins = .true.
596  DO ms = 1, mesh%mes
597  IF (.NOT.virgins(ms)) cycle
598  IF (minval(abs(mesh%sides(ms)-list_of_interfaces))/=0) cycle !==ms not on a cut
599  i_loc = mesh%jjs(:,ms)
600  DO msop = 1, mesh%mes
601  IF (msop==ms .OR. .NOT.virgins(msop)) cycle
602  IF (minval(abs(mesh%sides(msop)-list_of_interfaces))/=0) cycle !==msop not on a cut
603  DO ns = 1, nws
604  test = .false.
605  DO nsop = 1, nws
606  iop = mesh%jjs(nsop,msop)
607  IF (maxval(abs(mesh%rr(:,i_loc(ns))-mesh%rr(:,iop))).LT.epsilon) THEN
608  test = .true.
609  EXIT
610  END IF
611  END DO
612  IF (.NOT.test) THEN
613  EXIT !==This msop does not coincide with ms
614  END IF
615  END DO
616  IF (test) EXIT
617  END DO
618  IF (.NOT.test) THEN
619  CALL error_petsc(.NOT.'BUG in part_mesh_mhd, test ')
620  END IF
621  DO n = 1, 3
622  IF (neigh_new(n,mesh%neighs(msop))==0) THEN
623  neigh_new(n,mesh%neighs(msop)) = mesh%neighs(ms)
624  END IF
625  IF (neigh_new(n,mesh%neighs(ms))==0) THEN
626  neigh_new(n,mesh%neighs(ms)) = mesh%neighs(msop)
627  END IF
628  END DO
629  virgins(ms) = .false.
630  virgins(msop) = .false.
631  END DO
632  END IF
633 
634  !=================TEST PERIODIC==================
635  IF (PRESENT(my_periodic)) THEN
636  IF (my_periodic%nb_periodic_pairs/=0) THEN
637  DO ms = 1, mesh%mes
638  m = mesh%neighs(ms)
639  IF (minval(abs(mesh%sides(ms)-my_periodic%list_periodic(1,:))) == 0) THEN
640  jm_loc = minloc(abs(mesh%sides(ms)-my_periodic%list_periodic(1,:)))
641  s2 = my_periodic%list_periodic(2,jm_loc(1))
642  test = .false.
643  DO msop = 1, mesh%mes
644  IF (mesh%sides(msop) /= s2) cycle
645 
646  err = 0.d0
647  DO ns = 1, SIZE(my_periodic%vect_e,1)
648  err = err + abs(sum(mesh%rr(ns,mesh%jjs(:,ms))-mesh%rr(ns,mesh%jjs(:,msop)) &
649  +my_periodic%vect_e(ns,jm_loc(1))))
650  END DO
651 
652  IF (err .LE. epsilon) THEN
653  test = .true.
654  EXIT
655  END IF
656  END DO
657  IF (.NOT.test) THEN
658  CALL error_petsc('BUG, mop non trouve')
659  END IF
660  mop = mesh%neighs(msop)
661  DO n = 1, 3
662  IF (neigh_new(n,m) == 0) THEN
663  neigh_new(n,m) = mop
664  END IF
665  IF (neigh_new(n,mop) == 0) THEN
666  neigh_new(n,mop) = m
667  END IF
668  END DO
669  END IF
670  END DO
671  END IF
672  END IF
673  !================================================
674 
675 
676  ! Create the connectivity array based on neigh
677  nb_neigh = SIZE(mesh%neigh,1)
678  xind(1) = 1
679  DO m = 1, mesh%me
680  nb = 0
681  DO n = 1, nb_neigh
682  IF (neigh_new(n,m)==0) cycle
683  nb = nb + 1
684  END DO
685  xind(m+1) = xind(m) + nb
686  END DO
687  ALLOCATE(xadj(xind(mesh%me+1)-1))
688  p = 0
689  DO m = 1, mesh%me
690  DO n = 1, nb_neigh
691  IF (neigh_new(n,m)==0) cycle
692  p = p + 1
693  xadj(p) = neigh_new(n,m)
694  END DO
695  END DO
696  IF (p/=xind(mesh%me+1)-1) THEN
697  CALL error_petsc('BUG, p/=xind(mesh%me+1)-1')
698  END IF
699  ! End create the connectivity array based on neigh
700 
701  ALLOCATE(adjwgt(SIZE(xadj)))
702  opts = 0
703  !vwgt = 1
704  adjwgt = 1
705  numflag = 1 ! Fortran numbering of processors
706  wgtflag = 2
707  ALLOCATE(tpwgts(nb_proc))
708  tpwgts=one_k4/nb_proc
709  CALL metis_setdefaultoptions(metis_opt)
710  metis_opt(metis_option_numbering)=1
711  ubvec=1.01
712  CALL metis_partgraphrecursive(mesh%me, 1, xind,xadj,vwgt, vwgt, adjwgt, nb_proc,tpwgts , ubvec, metis_opt, edge, part)
713 !!$ CALL METIS_PartGraphRecursive(mesh%me,xind,xadj,vwgt,adjwgt,wgtflag, numflag, nb_proc, opts, edge, part)
714 
715  ! Create parts and modify part
716  !==Search on the boundary whether ms is on a cut.
717  nws = SIZE( mesh%jjs,1)
718  IF (SIZE(list_of_interfaces)/=0) THEN
719  virgins = .true.
720  DO ms = 1, mesh%mes
721  IF (.NOT.virgins(ms)) cycle
722  IF (minval(abs(mesh%sides(ms)-list_of_interfaces))/=0) cycle !==ms not on a cut
723  i_loc = mesh%jjs(:,ms)
724  DO msop = 1, mesh%mes
725  IF (msop==ms .OR. .NOT.virgins(msop)) cycle
726  IF (minval(abs(mesh%sides(msop)-list_of_interfaces))/=0) cycle !==msop not on a cut
727  DO ns = 1, nws
728  test = .false.
729  DO nsop = 1, nws
730  iop = mesh%jjs(nsop,msop)
731  IF (maxval(abs(mesh%rr(:,i_loc(ns))-mesh%rr(:,iop))).LT.epsilon) THEN
732  test = .true.
733  EXIT
734  END IF
735  END DO
736  IF (.NOT.test) THEN
737  EXIT !==This msop does not coincide with ms
738  END IF
739  END DO
740  IF (test) EXIT
741  END DO
742  IF (.NOT.test) THEN
743  CALL error_petsc(.NOT.'BUG in create_local_mesh, test ')
744  END IF
745  IF (part(mesh%neighs(ms)) == part(mesh%neighs(msop))) cycle !==ms is an internal cut
746  proc = min(part(mesh%neighs(ms)),part(mesh%neighs(msop)))
747  part(mesh%neighs(ms)) = proc !make sure interface are internal
748  part(mesh%neighs(msop)) = proc !make sure interface are internal
749  virgins(ms) = .false.
750  virgins(msop) = .false.
751  END DO
752  END IF
753  ! End create parts and modify part
754 
755  !=================TEST PERIODIC==================
756  IF (PRESENT(my_periodic)) THEN
757  IF (my_periodic%nb_periodic_pairs/=0) THEN
758 
759 
760  DO j = 1, mesh%np
761  per_pts(j,1) = j
762  END DO
763  per_pts(:,2:3) = 0
764  DO ms = 1, mesh%mes
765  m = mesh%neighs(ms)
766  IF ((minval(abs(mesh%sides(ms)-my_periodic%list_periodic(1,:))) /=0) .AND. &
767  (minval(abs(mesh%sides(ms)-my_periodic%list_periodic(2,:))) /=0) ) cycle
768  DO ns = 1, SIZE(mesh%jjs,1)
769  j = mesh%jjs(ns,ms)
770  per_pts(j,2) = m
771  DO msop = 1, mesh%mes
772  IF (minval(abs(mesh%sides(msop)-my_periodic%list_periodic(:,:))) /=0 ) cycle
773  IF (msop == ms) cycle
774  DO nsop = 1, SIZE(mesh%jjs,1)
775  IF (mesh%jjs(nsop,msop)==j) THEN
776  per_pts(j,3) = mesh%neighs(msop)
777  END IF
778  END DO
779  END DO
780  END DO
781  END DO
782 !!$ WRITE(*,*) per_pts(:,2)
783  CALL reassign_per_pts(mesh, part, per_pts)
784 
785  DO ms = 1, mesh%mes
786  m = mesh%neighs(ms)
787  IF (minval(abs(mesh%sides(ms)-my_periodic%list_periodic(1,:))) /= 0) cycle
788  jm_loc = minloc(abs(mesh%sides(ms)-my_periodic%list_periodic(1,:)))
789  s2 = my_periodic%list_periodic(2,jm_loc(1))
790  test = .false.
791  DO msop = 1, mesh%mes
792  IF (mesh%sides(msop) /= s2) cycle
793 
794  err = 0.d0
795  DO ns = 1, SIZE(my_periodic%vect_e,1)
796  err = err + abs(sum(mesh%rr(ns,mesh%jjs(:,ms))-mesh%rr(ns,mesh%jjs(:,msop)) &
797  +my_periodic%vect_e(ns,jm_loc(1))))
798  END DO
799 
800  IF (err .LE. epsilon) THEN
801  test = .true.
802  EXIT
803  END IF
804  END DO
805  IF (.NOT.test) THEN
806  CALL error_petsc('BUG, mop non trouve')
807  END IF
808  IF (part(mesh%neighs(ms)) /= part(mesh%neighs(msop))) THEN !==ms is an internal cut
809  proc = min(part(mesh%neighs(ms)),part(mesh%neighs(msop)))
810  part(mesh%neighs(ms)) = proc !make sure interface are internal
811  part(mesh%neighs(msop)) = proc !make sure interface are internal
812  END IF
813  END DO
814 
815  END IF
816  END IF
817 
818 !!$ WARNING, FL 1/2/13 : TO BE ADDED IN NEEDED
819  !================================================
820  CALL mpi_comm_rank(mpi_comm_world,rank,ierr)
821  IF (rank==0) THEN
822  CALL plot_const_p1_label(mesh%jj, mesh%rr, 1.d0*part, 'dd.plt')
823  END IF
824  !================================================
825 !!$ WARNING, FL 1/2/13 : TO BE ADDED IN NEEDED
826 
827  DEALLOCATE(xadj,adjwgt)
828  DEALLOCATE(tpwgts)
829 
830  END SUBROUTINE part_mesh_mhd
831 
832 
833  SUBROUTINE part_mesh_mhd_bis(nb_proc,list_u, list_h_in, list_phi ,mesh,list_of_interfaces,part,my_periodic)
834  USE def_type_mesh
835  USE my_util
836  USE sub_plot
837  USE periodic
838  IMPLICIT NONE
839  TYPE(mesh_type) :: mesh
840  INTEGER, DIMENSION(mesh%me) :: part
841  INTEGER, DIMENSION(:) :: list_of_interfaces
842  INTEGER, DIMENSION(:) :: list_u, list_h_in, list_phi
843  TYPE(periodic_data), OPTIONAL :: my_periodic
844 
845  LOGICAL, DIMENSION(mesh%mes) :: virgins
846  INTEGER, DIMENSION(3,mesh%me) :: neigh_new
847  INTEGER, DIMENSION(5) :: opts
848  INTEGER, DIMENSION(SIZE(mesh%jjs,1)) :: i_loc
849  INTEGER, DIMENSION(:), ALLOCATABLE :: xadj_u, xadj_h, xadj_phi, list_h
850  INTEGER, DIMENSION(:), ALLOCATABLE :: xind_u, xind_h, xind_phi
851  INTEGER, DIMENSION(:), ALLOCATABLE :: vwgt, adjwgt
852  INTEGER, DIMENSION(:), ALLOCATABLE :: u2glob, h2glob, phi2glob
853  INTEGER, DIMENSION(:), ALLOCATABLE :: part_u, part_h, part_phi
854  INTEGER, DIMENSION(1) :: jm_loc
855 !!$ INTEGER, DIMENSION(1) :: jloc
856  INTEGER, DIMENSION(mesh%np,3) :: per_pts
857  INTEGER, DIMENSION(mesh%me) :: glob2loc
858  INTEGER, DIMENSION(mesh%np) :: indicator
859  INTEGER, DIMENSION(3) :: j_loc
860  INTEGER :: nb_neigh, edge, m, ms, n, nb, numflag, p, wgtflag, j, &
861  ns, nws, msop, nsop, proc, iop, mop, s2, k, me_u, me_h, me_phi, idm
862  REAL(KIND=8) :: err
863  LOGICAL :: test
864  !===(JLG) Feb 20, 2019. Petsc developpers decide to use REAL(KIND=4) to interface with metis
865  !REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: tpwgts
866  !REAL(KIND=8), DIMENSION(1) :: ubvec
867  REAL(KIND=4), DIMENSION(:), ALLOCATABLE :: tpwgts
868  REAL(KIND=4), DIMENSION(1) :: ubvec
869  REAL(KIND=4) :: one_K4=1.0
870  !===(JLG)Feb 20, 2019.
871  INTEGER, DIMENSION(METIS_NOPTIONS) :: metis_opt
872  petscmpiint :: nb_proc
873 !!$ WARNING, FL 1/2/13 : TO BE ADDED IN NEEDED
874  petscerrorcode :: ierr
875  petscmpiint :: rank
876 !!$ WARNING, FL 1/2/13 : TO BE ADDED IN NEEDED
877 
878  IF (nb_proc==1) THEN
879  part = 1
880  RETURN
881  END IF
882  glob2loc = 0
883  !Create list_h = list_h_in \ list_ns
884  nb = 0
885  DO j = 1, SIZE(list_h_in)
886  IF (minval(abs(list_h_in(j)-list_u))==0) cycle
887  nb = nb + 1
888  END DO
889  ALLOCATE(list_h(nb))
890  nb = 0
891  DO j = 1, SIZE(list_h_in)
892  IF (minval(abs(list_h_in(j)-list_u))==0) cycle
893  nb = nb + 1
894  list_h(nb) = list_h_in(j)
895  END DO
896 
897  !Done with list_h
898 
899  nws = SIZE( mesh%jjs,1)
900  neigh_new = mesh%neigh
901  IF (SIZE(list_of_interfaces)/=0) THEN
902  virgins = .true.
903  DO ms = 1, mesh%mes
904  IF (.NOT.virgins(ms)) cycle
905  IF (minval(abs(mesh%sides(ms)-list_of_interfaces))/=0) cycle !==ms not on a cut
906  i_loc = mesh%jjs(:,ms)
907  DO msop = 1, mesh%mes
908  IF (msop==ms .OR. .NOT.virgins(msop)) cycle
909  IF (minval(abs(mesh%sides(msop)-list_of_interfaces))/=0) cycle !==msop not on a cut
910  DO ns = 1, nws
911  test = .false.
912  DO nsop = 1, nws
913  iop = mesh%jjs(nsop,msop)
914  IF (maxval(abs(mesh%rr(:,i_loc(ns))-mesh%rr(:,iop))).LT.epsilon) THEN
915  test = .true.
916  EXIT
917  END IF
918  END DO
919  IF (.NOT.test) THEN
920  EXIT !==This msop does not coincide with ms
921  END IF
922  END DO
923  IF (test) EXIT
924  END DO
925  IF (.NOT.test) THEN
926  CALL error_petsc(.NOT.'BUG in part_mesh_mhd, test ')
927  END IF
928  DO n = 1, 3
929  IF (neigh_new(n,mesh%neighs(msop))==0) THEN
930  neigh_new(n,mesh%neighs(msop)) = mesh%neighs(ms)
931  END IF
932  IF (neigh_new(n,mesh%neighs(ms))==0) THEN
933  neigh_new(n,mesh%neighs(ms)) = mesh%neighs(msop)
934  END IF
935  END DO
936  virgins(ms) = .false.
937  virgins(msop) = .false.
938  END DO
939  END IF
940 
941  !=================TEST PERIODIC==================
942 
943  IF (PRESENT(my_periodic)) THEN
944  IF (my_periodic%nb_periodic_pairs/=0) THEN
945  DO ms = 1, mesh%mes
946  m = mesh%neighs(ms)
947  IF (minval(abs(mesh%sides(ms)-my_periodic%list_periodic(1,:))) == 0) THEN
948  jm_loc = minloc(abs(mesh%sides(ms)-my_periodic%list_periodic(1,:)))
949  s2 = my_periodic%list_periodic(2,jm_loc(1))
950  test = .false.
951  DO msop = 1, mesh%mes
952  IF (mesh%sides(msop) /= s2) cycle
953 
954  err = 0.d0
955  DO ns = 1, SIZE(my_periodic%vect_e,1)
956  err = err + abs(sum(mesh%rr(ns,mesh%jjs(:,ms))-mesh%rr(ns,mesh%jjs(:,msop)) &
957  +my_periodic%vect_e(ns,jm_loc(1))))
958  END DO
959 
960  IF (err .LE. epsilon) THEN
961  test = .true.
962  EXIT
963  END IF
964  END DO
965  IF (.NOT.test) THEN
966  CALL error_petsc('BUG, mop non trouve')
967  END IF
968  mop = mesh%neighs(msop)
969  DO n = 1, 3
970  IF (neigh_new(n,m) == 0) THEN
971  neigh_new(n,m) = mop
972  END IF
973  IF (neigh_new(n,mop) == 0) THEN
974  neigh_new(n,mop) = m
975  END IF
976  END DO
977  END IF
978  END DO
979  END IF
980  END IF
981 
982  !================================================
983  me_u = 0
984  me_h = 0
985  me_phi = 0
986  DO m = 1, mesh%me
987  idm = mesh%i_d(m)
988  IF (minval(abs(idm-list_u))==0) THEN
989  me_u = me_u + 1
990  ELSE IF (minval(abs(idm-list_h))==0) THEN
991  me_h = me_h + 1
992  ELSE IF (minval(abs(idm-list_phi))==0) THEN
993  me_phi = me_phi + 1
994  ELSE
995  CALL error_petsc('BUG in part_mesh_mhd_bis : element not in the mesh')
996  END IF
997  END DO
998  ALLOCATE(u2glob(me_u), h2glob(me_h),phi2glob(me_phi))
999  me_u = 0
1000  me_h = 0
1001  me_phi = 0
1002  DO m = 1, mesh%me
1003  idm = mesh%i_d(m)
1004  IF (minval(abs(idm-list_u))==0) THEN
1005  me_u = me_u + 1
1006  u2glob(me_u) = m
1007  glob2loc(m) = me_u
1008  ELSE IF (minval(abs(idm-list_h))==0) THEN
1009  me_h = me_h + 1
1010  h2glob(me_h) = m
1011  glob2loc(m) = me_h
1012  ELSE IF (minval(abs(idm-list_phi))==0) THEN
1013  me_phi = me_phi + 1
1014  phi2glob(me_phi) = m
1015  glob2loc(m) = me_phi
1016  ELSE
1017  CALL error_petsc('BUG in part_mesh_mhd_bis : element not in the mesh')
1018  END IF
1019  END DO
1020  ! Create the connectivity arrays based on neigh
1021  nb_neigh = SIZE(mesh%neigh,1)
1022  ALLOCATE(xind_u(me_u+1), xind_h(me_h+1), xind_phi(me_phi+1))
1023  xind_u(1) = 1
1024  DO k = 1, me_u
1025  m = u2glob(k)
1026  nb = 0
1027  DO n = 1, nb_neigh
1028  mop = neigh_new(n,m)
1029  IF (mop==0) cycle
1030  IF (minval(abs(mesh%i_d(mop)-list_u))/=0) cycle
1031  nb = nb + 1
1032  END DO
1033  xind_u(k+1) = xind_u(k) + nb
1034  END DO
1035  xind_h(1) = 1
1036  DO k = 1, me_h
1037  m = h2glob(k)
1038  nb = 0
1039  DO n = 1, nb_neigh
1040  mop = neigh_new(n,m)
1041  IF (mop==0) cycle
1042  IF (minval(abs(mesh%i_d(mop)-list_h))/=0) cycle
1043  nb = nb + 1
1044  END DO
1045  xind_h(k+1) = xind_h(k) + nb
1046  END DO
1047  xind_phi(1) = 1
1048  DO k = 1, me_phi
1049  m = phi2glob(k)
1050  nb = 0
1051  DO n = 1, nb_neigh
1052  mop = neigh_new(n,m)
1053  IF (mop==0) cycle
1054  IF (minval(abs(mesh%i_d(mop)-list_phi))/=0) cycle
1055  nb = nb + 1
1056  END DO
1057  xind_phi(k+1) = xind_phi(k) + nb
1058  END DO
1059 
1060  ALLOCATE(xadj_u(xind_u(me_u+1)-1))
1061  ALLOCATE(xadj_h(xind_h(me_h+1)-1))
1062  ALLOCATE(xadj_phi(xind_phi(me_phi+1)-1))
1063  p = 0
1064  DO k = 1, me_u
1065  m = u2glob(k)
1066  DO n = 1, nb_neigh
1067  mop = neigh_new(n,m)
1068  IF (mop==0) cycle
1069  IF (minval(abs(mesh%i_d(mop)-list_u))/=0) cycle
1070 !!$ jloc = MINLOC(ABS(u2glob-mop))
1071  p = p + 1
1072 !!$ xadj_u(p) = jloc(1)
1073  xadj_u(p) = glob2loc(mop)
1074  END DO
1075  END DO
1076  IF (p/=xind_u(me_u+1)-1) THEN
1077  CALL error_petsc('BUG, p/=xind_u(me_u+1)-1')
1078  END IF
1079  p = 0
1080  DO k = 1, me_h
1081  m = h2glob(k)
1082  DO n = 1, nb_neigh
1083  mop = neigh_new(n,m)
1084  IF (mop==0) cycle
1085  IF (minval(abs(mesh%i_d(mop)-list_h))/=0) cycle
1086 !!$ jloc = MINLOC(ABS(h2glob-mop))
1087  p = p + 1
1088 !!$ xadj_h(p) = jloc(1)
1089  xadj_h(p) = glob2loc(mop)
1090  END DO
1091  END DO
1092  IF (p/=xind_h(me_h+1)-1) THEN
1093  CALL error_petsc('BUG, p/=xind_h(me_h+1)-1')
1094  END IF
1095  p = 0
1096  DO k = 1, me_phi
1097  m = phi2glob(k)
1098  DO n = 1, nb_neigh
1099  mop = neigh_new(n,m)
1100  IF (mop==0) cycle
1101  IF (minval(abs(mesh%i_d(mop)-list_phi))/=0) cycle
1102 !!$ jloc = MINLOC(ABS(phi2glob-mop))
1103  p = p + 1
1104 !!$ xadj_phi(p) = jloc(1)
1105  xadj_phi(p) = glob2loc(mop)
1106  END DO
1107  END DO
1108  IF (p/=xind_phi(me_phi+1)-1) THEN
1109  CALL error_petsc('BUG, p/=xind_phi(me_phi+1)-1')
1110  END IF
1111 
1112  opts = 0
1113  numflag = 1
1114  wgtflag = 2
1115  ! Create partitions
1116  ALLOCATE(tpwgts(nb_proc))
1117  tpwgts=one_k4/nb_proc
1118  CALL metis_setdefaultoptions(metis_opt)
1119  metis_opt(metis_option_numbering)=1
1120  ubvec=1.01
1121  IF (me_u /= 0) THEN
1122  ALLOCATE(vwgt(me_u), adjwgt(SIZE(xadj_u)), part_u(me_u))
1123  vwgt = 1
1124  adjwgt = 1
1125  CALL metis_partgraphrecursive(me_u, 1, xind_u,xadj_u,vwgt, vwgt, adjwgt, nb_proc,tpwgts , ubvec, metis_opt, edge, part_u)
1126 !!$ CALL METIS_PartGraphRecursive(me_u,xind_u,xadj_u,vwgt,adjwgt,wgtflag, numflag, nb_proc, opts, edge, part_u)
1127  END IF
1128  IF (me_h /= 0) THEN
1129  IF (ALLOCATED(vwgt)) THEN
1130  DEALLOCATE(vwgt, adjwgt)
1131  END IF
1132  ALLOCATE(vwgt(me_h), adjwgt(SIZE(xadj_h)), part_h(me_h))
1133  vwgt = 1
1134  adjwgt = 1
1135  CALL metis_partgraphrecursive(me_h, 1, xind_h,xadj_h,vwgt, vwgt, adjwgt, nb_proc,tpwgts , ubvec, metis_opt, edge, part_h)
1136 !!$ CALL METIS_PartGraphRecursive(me_h,xind_h,xadj_h,vwgt,adjwgt,wgtflag, numflag, nb_proc, opts, edge, part_h)
1137  END IF
1138  IF (me_phi /= 0) THEN
1139  IF (ALLOCATED(vwgt)) THEN
1140  DEALLOCATE(vwgt, adjwgt)
1141  END IF
1142  ALLOCATE(vwgt(me_phi), adjwgt(SIZE(xadj_phi)), part_phi(me_phi))
1143  vwgt = 1
1144  adjwgt = 1
1145  CALL metis_partgraphrecursive(me_phi, 1, xind_phi,xadj_phi,vwgt, vwgt, adjwgt, nb_proc,tpwgts, &
1146  ubvec, metis_opt, edge, part_phi)
1147 !!$ CALL METIS_PartGraphRecursive(me_phi,xind_phi,xadj_phi,vwgt,adjwgt,wgtflag, numflag, nb_proc, opts, edge, part_phi)
1148  END IF
1149 !!$!TESTTT
1150 !!$ do m = 1, mesh%me
1151 !!$ if (part_phi(m) ==1) THEN
1152 !!$ part_phi(m) =2
1153 !!$ ELSE if (part_phi(m) ==2) then
1154 !!$ part_phi(m) =1
1155 !!$ ELSE if (part_phi(m) ==3) then
1156 !!$ part_phi(m) =4
1157 !!$ ELSE if (part_phi(m) ==4) then
1158 !!$ part_phi(m) =3
1159 !!$ END IF
1160 !!$ end do
1161 !!$!TESTTT
1162  ! Create global partition 'part'
1163  part = -1
1164  IF (me_u/=0) THEN
1165  part(u2glob(:)) = part_u
1166  END IF
1167  IF (me_h/=0) THEN
1168  part(h2glob(:)) = part_h
1169  END IF
1170  IF (me_phi/=0) THEN
1171  part(phi2glob(:)) = part_phi
1172  END IF
1173  IF (minval(part)==-1) THEN
1174  CALL error_petsc('BUG in part_mesh_mhd_bis, MINVAL(part) == -1')
1175  END IF
1176  ! End Create global partition
1177 
1178  ! Create parts and modify part
1179  !==Search on the boundary whether ms is on a cut.
1180  IF (SIZE(mesh%jj,1)/=3) THEN
1181  write(*,*) 'SIZE(mesh%jj,1)', SIZE(mesh%jj,1)
1182  CALL error_petsc('BUG in part_mesh_mhd_bis, SIZE(mesh%jj,1)/=3')
1183  END IF
1184  indicator = -1
1185  nws = SIZE( mesh%jjs,1)
1186  IF (SIZE(list_of_interfaces)/=0) THEN
1187  virgins = .true.
1188  DO ms = 1, mesh%mes
1189  IF (.NOT.virgins(ms)) cycle
1190  IF (minval(abs(mesh%sides(ms)-list_of_interfaces))/=0) cycle !==ms not on a cut
1191  i_loc = mesh%jjs(:,ms)
1192  DO msop = 1, mesh%mes
1193  IF (msop==ms .OR. .NOT.virgins(msop)) cycle
1194  IF (minval(abs(mesh%sides(msop)-list_of_interfaces))/=0) cycle !==msop not on a cut
1195  DO ns = 1, nws
1196  test = .false.
1197  DO nsop = 1, nws
1198  iop = mesh%jjs(nsop,msop)
1199  IF (maxval(abs(mesh%rr(:,i_loc(ns))-mesh%rr(:,iop))).LT.epsilon) THEN
1200  test = .true.
1201  EXIT
1202  END IF
1203  END DO
1204  IF (.NOT.test) THEN
1205  EXIT !==This msop does not coincide with ms
1206  END IF
1207  END DO
1208  IF (test) EXIT
1209  END DO
1210  IF (.NOT.test) THEN
1211  CALL error_petsc(.NOT.'BUG in create_local_mesh, test ')
1212  END IF
1213  IF (part(mesh%neighs(ms)) == part(mesh%neighs(msop))) cycle !==ms is an internal cut
1214  proc = min(part(mesh%neighs(ms)),part(mesh%neighs(msop)))
1215  part(mesh%neighs(ms)) = proc !make sure interface are internal
1216  part(mesh%neighs(msop)) = proc !make sure interface are internal
1217  virgins(ms) = .false.
1218  virgins(msop) = .false.
1219  indicator(mesh%jjs(:,ms)) = proc
1220  indicator(mesh%jjs(:,msop)) = proc
1221  END DO
1222  END IF
1223  !===Fix the partition so that all the cells having one vertex on an
1224  !===interface belong to the same processor as those sharing this vertices and
1225  !===having two vertices on the interface (JLG + DCL July 22 2015)
1226  DO m = 1, mesh%me
1227  j_loc = mesh%jj(:,m)
1228  n = maxval(indicator(j_loc))
1229  IF (n == -1) cycle
1230  IF (indicator(j_loc(1))*indicator(j_loc(2))*indicator(j_loc(3))<0) cycle
1231  part(m) = n
1232  END DO
1233  ! End create parts and modify part
1234 
1235  !=================TEST PERIODIC==================
1236  IF (PRESENT(my_periodic)) THEN
1237  IF (my_periodic%nb_periodic_pairs/=0) THEN
1238 
1239 
1240  DO j = 1, mesh%np
1241  per_pts(j,1) = j
1242  END DO
1243  per_pts(:,2:3) = 0
1244  DO ms = 1, mesh%mes
1245  m = mesh%neighs(ms)
1246  IF ((minval(abs(mesh%sides(ms)-my_periodic%list_periodic(1,:))) /=0) .AND. &
1247  (minval(abs(mesh%sides(ms)-my_periodic%list_periodic(2,:))) /=0) ) cycle
1248  DO ns = 1, SIZE(mesh%jjs,1)
1249  j = mesh%jjs(ns,ms)
1250  per_pts(j,2) = m
1251  DO msop = 1, mesh%mes
1252  IF (minval(abs(mesh%sides(msop)-my_periodic%list_periodic(:,:))) /=0 ) cycle
1253  IF (msop == ms) cycle
1254  DO nsop = 1, SIZE(mesh%jjs,1)
1255  IF (mesh%jjs(nsop,msop)==j) THEN
1256  per_pts(j,3) = mesh%neighs(msop)
1257  END IF
1258  END DO
1259  END DO
1260  END DO
1261  END DO
1262 !!$ WRITE(*,*) per_pts(:,2)
1263  CALL reassign_per_pts(mesh, part, per_pts)
1264 
1265  DO ms = 1, mesh%mes
1266  m = mesh%neighs(ms)
1267  IF (minval(abs(mesh%sides(ms)-my_periodic%list_periodic(1,:))) /= 0) cycle
1268  jm_loc = minloc(abs(mesh%sides(ms)-my_periodic%list_periodic(1,:)))
1269  s2 = my_periodic%list_periodic(2,jm_loc(1))
1270  test = .false.
1271  DO msop = 1, mesh%mes
1272  IF (mesh%sides(msop) /= s2) cycle
1273 
1274  err = 0.d0
1275  DO ns = 1, SIZE(my_periodic%vect_e,1)
1276  err = err + abs(sum(mesh%rr(ns,mesh%jjs(:,ms))-mesh%rr(ns,mesh%jjs(:,msop)) &
1277  +my_periodic%vect_e(ns,jm_loc(1))))
1278  END DO
1279 
1280  IF (err .LE. epsilon) THEN
1281  test = .true.
1282  EXIT
1283  END IF
1284  END DO
1285  IF (.NOT.test) THEN
1286  CALL error_petsc('BUG, mop non trouve')
1287  END IF
1288  IF (part(mesh%neighs(ms)) /= part(mesh%neighs(msop))) THEN !==ms is an internal cut
1289  proc = min(part(mesh%neighs(ms)),part(mesh%neighs(msop)))
1290  part(mesh%neighs(ms)) = proc !make sure interface are internal
1291  part(mesh%neighs(msop)) = proc !make sure interface are internal
1292  END IF
1293  END DO
1294 
1295  END IF
1296  END IF
1297 
1298 !!$ WARNING, FL 1/2/13 : TO BE ADDED IN NEEDED
1299  !================================================
1300  CALL mpi_comm_rank(mpi_comm_world,rank,ierr)
1301  IF (rank==0) THEN
1302  CALL plot_const_p1_label(mesh%jj, mesh%rr, 1.d0*part, 'dd.plt')
1303  END IF
1304  !================================================
1305 !!$ WARNING, FL 1/2/13 : TO BE ADDED IN NEEDED
1306 
1307  DEALLOCATE(vwgt,adjwgt)
1308  IF (ALLOCATED(xadj_u)) DEALLOCATE(xadj_u)
1309  IF (ALLOCATED(xadj_h)) DEALLOCATE(xadj_h)
1310  IF (ALLOCATED(xadj_phi)) DEALLOCATE(xadj_phi)
1311  IF (ALLOCATED(list_h)) DEALLOCATE(list_h)
1312  IF (ALLOCATED(xind_u)) DEALLOCATE(xind_u)
1313  IF (ALLOCATED(xind_h)) DEALLOCATE(xind_h)
1314  IF (ALLOCATED(xind_phi)) DEALLOCATE(xind_phi)
1315  IF (ALLOCATED(u2glob)) DEALLOCATE(u2glob)
1316  IF (ALLOCATED(h2glob)) DEALLOCATE(h2glob)
1317  IF (ALLOCATED(phi2glob)) DEALLOCATE(phi2glob)
1318  IF (ALLOCATED(part_u)) DEALLOCATE(part_u)
1319  IF (ALLOCATED(part_h)) DEALLOCATE(part_h)
1320  IF (ALLOCATED(part_phi)) DEALLOCATE(part_phi)
1321 
1322  DEALLOCATE(tpwgts)
1323 
1324  END SUBROUTINE part_mesh_mhd_bis
1325 
1326  SUBROUTINE extract_mesh(communicator,nb_proc,mesh_glob,part,list_dom,mesh,mesh_loc)
1327  USE def_type_mesh
1328  USE my_util
1329  IMPLICIT NONE
1330  TYPE(mesh_type) :: mesh_glob, mesh, mesh_loc
1331  INTEGER, DIMENSION(:) :: part, list_dom
1332  INTEGER, DIMENSION(mesh_glob%me) :: bat
1333  INTEGER, DIMENSION(mesh_glob%np) :: i_old_to_new
1334  INTEGER, DIMENSION(mesh_glob%mes) :: parts
1335  INTEGER, DIMENSION(nb_proc) :: nblmt_per_proc, start, displ
1336  INTEGER, DIMENSION(2) :: np_loc, me_loc, mes_loc
1337  INTEGER, DIMENSION(:), ALLOCATABLE :: list_m, tab, tabs
1338  INTEGER :: nb_proc, ms, i, index, m, mop, n
1339  petscerrorcode :: ierr
1340  petscmpiint :: rank
1341  mpi_comm :: communicator
1342  CALL mpi_comm_rank(communicator,rank,ierr)
1343 
1344  ! Create parts
1345  parts = part(mesh_glob%neighs)
1346  ! End create parts
1347 
1348  ! Create list_m
1349  i = 0
1350  DO m = 1, mesh_glob%me
1351  IF (minval(abs(list_dom-mesh_glob%i_d(m)))/=0) cycle
1352  i = i + 1
1353  END DO
1354  mesh%me = i
1355  ALLOCATE (list_m(mesh%me))
1356  i = 0
1357  DO m = 1, mesh_glob%me
1358  IF (minval(abs(list_dom-mesh_glob%i_d(m)))/=0) cycle
1359  i = i + 1
1360  list_m(i) = m
1361  END DO
1362  !End create list_m
1363 
1364  ! Count elements on processors
1365  nblmt_per_proc = 0
1366  DO i = 1, mesh%me
1367  m = list_m(i)
1368  nblmt_per_proc(part(m)) = nblmt_per_proc(part(m)) + 1
1369  END DO
1370  start(1) = 0
1371  DO n = 2, nb_proc
1372  start(n) = start(n-1) + nblmt_per_proc(n-1)
1373  END DO
1374  me_loc(1) = start(rank+1) + 1
1375  me_loc(2) = start(rank+1) + nblmt_per_proc(rank+1)
1376  displ = start
1377  ! End count elements on processors
1378 
1379  ! Re-order elements
1380  ALLOCATE(tab(mesh%me))
1381  bat = 0
1382  DO i = 1, mesh%me
1383  m = list_m(i)
1384  start(part(m)) = start(part(m)) + 1
1385  tab(start(part(m))) = m
1386  bat(m) = start(part(m))
1387  END DO
1388  ! Re-order elements
1389 
1390  ! Create mesh%jj
1391  mesh%gauss%n_w = SIZE(mesh_glob%jj,1)
1392  ALLOCATE(mesh%jj(SIZE(mesh_glob%jj,1),mesh%me))
1393  i_old_to_new = 0
1394  index = 0
1395  DO m = 1, mesh%me
1396  DO n = 1, SIZE(mesh_glob%jj,1)
1397  i = mesh_glob%jj(n,tab(m))
1398  IF (i_old_to_new(i)/=0) THEN
1399  mesh%jj(n,m) = i_old_to_new(i)
1400  ELSE
1401  index = index + 1
1402  i_old_to_new(i) = index
1403  mesh%jj(n,m) = i_old_to_new(i)
1404  END IF
1405  END DO
1406  END DO
1407  ! End Create mesh%jj
1408 
1409  ! Create mesh%rr
1410  mesh%np = index
1411  ALLOCATE(mesh%rr(2,mesh%np))
1412  DO i = 1, mesh_glob%np
1413  IF (i_old_to_new(i)==0) cycle
1414  mesh%rr(:,i_old_to_new(i)) = mesh_glob%rr(:,i)
1415  END DO
1416  !End Create mesh%rr
1417 
1418  ! Create mesh%neigh
1419  ALLOCATE(mesh%neigh(3,mesh%me))
1420  DO m = 1, mesh%me
1421  DO n = 1, 3
1422  mop = mesh_glob%neigh(n,tab(m))
1423  IF (mop==0) THEN
1424  mesh%neigh(n,m) = 0
1425  ELSE
1426  mesh%neigh(n,m) = bat(mop)
1427  END IF
1428  END DO
1429  END DO
1430  ! End Create mesh%neigh
1431 
1432  ! Create mesh%i_d
1433  ALLOCATE(mesh%i_d(mesh%me))
1434  mesh%i_d = mesh_glob%i_d(tab)
1435  ! End mesh%i_d
1436 
1437  ! Create np_loc
1438  IF (displ(rank+1)/=0) THEN
1439  np_loc(1) = maxval(mesh%jj(:,1:displ(rank+1))) + 1
1440  ELSE
1441  np_loc(1) = 1
1442  END IF
1443  np_loc(2) = np_loc(1) - 1
1444  IF (me_loc(1).LE.me_loc(2)) THEN
1445  np_loc(2) = maxval(mesh%jj(:,me_loc(1):me_loc(2)))
1446  END IF
1447  IF (np_loc(2) .LT. np_loc(1)-1) THEN
1448  np_loc(2) = np_loc(1) - 1
1449  END IF
1450  ! End create np_loc
1451 
1452  ! Create mes_loc
1453  nblmt_per_proc=0
1454  DO ms = 1, mesh_glob%mes
1455  IF (minval(abs(list_dom-mesh_glob%i_d(mesh_glob%neighs(ms))))/=0) cycle
1456  n = parts(ms)
1457  nblmt_per_proc(n) = nblmt_per_proc(n) + 1
1458  END DO
1459  start(1) = 0
1460  DO n = 2, nb_proc
1461  start(n) = start(n-1) + nblmt_per_proc(n-1)
1462  END DO
1463  mes_loc(1) = start(rank+1) + 1
1464  mes_loc(2) = start(rank+1) + nblmt_per_proc(rank+1)
1465  mesh%mes = sum(nblmt_per_proc)
1466  ! End create mes_loc
1467 
1468  ! Create tabs and sbat
1469  ALLOCATE(tabs(mesh%mes))
1470  DO ms = 1, mesh_glob%mes
1471  IF (minval(abs(list_dom-mesh_glob%i_d(mesh_glob%neighs(ms))))/=0) cycle
1472  start(parts(ms)) = start(parts(ms)) + 1
1473  tabs(start(parts(ms))) = ms
1474  END DO
1475  ! End create tabs and sbat
1476 
1477  ! Create neighs
1478  ALLOCATE(mesh%neighs(mesh%mes))
1479  mesh%neighs = bat(mesh_glob%neighs(tabs))
1480  ! End create neighs
1481 
1482  ! Re-order sides
1483  ALLOCATE(mesh%sides(mesh%mes))
1484  mesh%sides = mesh_glob%sides(tabs)
1485  ! End re-order sides
1486 
1487  ! Re-order jjs
1488  mesh%gauss%n_ws = SIZE(mesh_glob%jjs,1)
1489  ALLOCATE(mesh%jjs(SIZE(mesh_glob%jjs,1),mesh%mes))
1490 
1491  DO n = 1, SIZE(mesh%jjs,1)
1492  mesh%jjs(n,:) = i_old_to_new(mesh_glob%jjs(n,tabs))
1493  END DO
1494  ! End re-order jjs
1495 
1496  !==We create the local mesh now
1497  mesh%edge_stab = .false.
1498  CALL create_local_mesh(mesh, mesh_loc, me_loc, mes_loc, np_loc)
1499 
1500  DEALLOCATE(list_m, tab, tabs)
1501 
1502  END SUBROUTINE extract_mesh
1503 
1504  SUBROUTINE create_local_mesh(mesh, mesh_loc, me_loc, mes_loc, np_loc)
1505  USE def_type_mesh
1506  USE my_util
1507  IMPLICIT NONE
1508  TYPE(mesh_type) :: mesh, mesh_loc
1509  INTEGER, DIMENSION(2), INTENT(IN) :: me_loc, mes_loc, np_loc
1510  INTEGER, DIMENSION(mesh%me) :: m_glob_to_loc
1511  INTEGER, DIMENSION(:), ALLOCATABLE :: m_loc_to_glob
1512  INTEGER, DIMENSION(mesh%np) :: glob_to_loc,loc_to_glob
1513  LOGICAL, DIMENSION(mesh%np) :: virgin
1514  INTEGER :: dim, nws, nw, m, ms, mop, ns, msup, minf, dof, &
1515  dom_me, nwc, dom_mes, dom_np, n, i
1516  LOGICAL :: test
1517 
1518  dim = SIZE(mesh%rr,1)
1519  nws = SIZE(mesh%jjs,1)
1520  nw = SIZE(mesh%jj,1)
1521  nwc = SIZE(mesh%neigh,1)
1522  !==Test if one proc only
1523  IF (me_loc(2) - me_loc(1) + 1==mesh%me) THEN
1524  mesh_loc%me = mesh%me
1525  mesh_loc%np = mesh%np
1526  mesh_loc%mes = mesh%mes
1527  mesh_loc%dom_me = mesh%me
1528  mesh_loc%dom_np = mesh%np
1529  mesh_loc%dom_mes = mesh%mes
1530  mesh_loc%gauss%n_w = nw
1531  ALLOCATE(mesh_loc%jj(nw,mesh%me))
1532  mesh_loc%jj = mesh%jj
1533  ALLOCATE(mesh_loc%loc_to_glob(mesh%np))
1534  DO n = 1, mesh%np
1535  mesh_loc%loc_to_glob(n) = n
1536  END DO
1537  ALLOCATE(mesh_loc%rr(dim,mesh%np))
1538  mesh_loc%rr=mesh%rr
1539  ALLOCATE(mesh_loc%neigh(nwc,mesh%me))
1540  mesh_loc%neigh = mesh%neigh
1541  ALLOCATE(mesh_loc%i_d(mesh%me))
1542  mesh_loc%i_d = mesh%i_d
1543  ALLOCATE(mesh_loc%neighs(mesh_loc%mes))
1544  mesh_loc%neighs = mesh%neighs
1545  ALLOCATE(mesh_loc%sides(mesh_loc%mes))
1546  mesh_loc%sides = mesh%sides
1547  mesh_loc%gauss%n_ws = nws
1548  ALLOCATE(mesh_loc%jjs(nws,mesh_loc%mes))
1549  mesh_loc%jjs = mesh%jjs
1550  RETURN
1551  END IF
1552  !==End test if one proc only
1553 
1554 
1555  !==Create the new mesh
1556  dom_me = me_loc(2) - me_loc(1) + 1
1557  dom_mes = mes_loc(2) - mes_loc(1) + 1
1558  dom_np = np_loc(2) - np_loc(1) + 1
1559  mesh_loc%me = dom_me
1560  mesh_loc%mes = dom_mes
1561  mesh_loc%dom_me = dom_me
1562  mesh_loc%dom_np = dom_np
1563  mesh_loc%dom_mes = dom_mes
1564 
1565  IF (dom_me==0) THEN
1566  ALLOCATE(mesh_loc%jj(0,0))
1567  ALLOCATE(mesh_loc%loc_to_glob(0))
1568  ALLOCATE(mesh_loc%rr(0,0))
1569  ALLOCATE(mesh_loc%neigh(0,0))
1570  ALLOCATE(mesh_loc%i_d(0))
1571  ALLOCATE(mesh_loc%neighs(0))
1572  ALLOCATE(mesh_loc%sides(0))
1573  ALLOCATE(mesh_loc%jjs(0,0))
1574  mesh_loc%gauss%n_w = 0
1575  mesh_loc%gauss%n_ws = 0
1576  RETURN
1577  ELSE IF (dom_me<0) THEN
1578  CALL error_petsc('BUG in create_local_mesh, dom_me<0 ')
1579  ELSE
1580  mesh_loc%gauss%n_w = nw
1581  mesh_loc%gauss%n_ws = nws
1582  END IF
1583 
1584  !==Re-order jj
1585  virgin = .true.
1586  dof = 0
1587  DO m = me_loc(1), me_loc(2)
1588  DO n = 1, nw
1589  i = mesh%jj(n,m)
1590  IF(.NOT.virgin(i) .OR. i.GE.np_loc(1)) cycle
1591  virgin(i) = .false.
1592  dof = dof + 1
1593  END DO
1594  END DO
1595  ALLOCATE(mesh_loc%jj(nw,mesh_loc%me))
1596 
1597  ALLOCATE(m_loc_to_glob(mesh_loc%me))
1598  m_glob_to_loc = 0
1599  virgin = .true.
1600  dof = dom_np
1601  DO m = me_loc(1), me_loc(2)
1602  DO n = 1, nw
1603  i = mesh%jj(n,m)
1604  IF(virgin(i)) THEN
1605  virgin(i) = .false.
1606  IF (i .LT. np_loc(1)) THEN
1607  dof = dof + 1
1608  glob_to_loc(i) = dof
1609  loc_to_glob(dof) = i
1610  ELSE
1611  glob_to_loc(i) = i-np_loc(1) + 1
1612  loc_to_glob(i-np_loc(1) + 1) = i
1613  END IF
1614  END IF
1615  END DO
1616  m_loc_to_glob(m-me_loc(1)+1) = m
1617  m_glob_to_loc(m) = m-me_loc(1)+1
1618  END DO
1619 
1620  DO n = 1, nw
1621  mesh_loc%jj(n,1:dom_me) = glob_to_loc(mesh%jj(n,me_loc(1):me_loc(2)))
1622  END DO
1623  !==End re-order jj
1624 
1625  !==Create mesh%loc_to_glob
1626  IF (maxval(mesh_loc%jj)/=dof) THEN
1627  CALL error_petsc('BUG in create_local_mesh, mesh_loc%jj)/=dof')
1628  END IF
1629  mesh_loc%np = dof
1630  ALLOCATE(mesh_loc%loc_to_glob(mesh_loc%np))
1631  mesh_loc%loc_to_glob = loc_to_glob(1:mesh_loc%np)
1632  !==End create mesh%loc_to_glob
1633 
1634  !==Re-order rr
1635  ALLOCATE(mesh_loc%rr(dim,mesh_loc%np))
1636  DO n = 1, mesh_loc%np
1637  mesh_loc%rr(:,n) = mesh%rr(:,mesh_loc%loc_to_glob(n))
1638  END DO
1639  !==End re-order rr
1640 
1641  !==Re-order neigh
1642  ALLOCATE(mesh_loc%neigh(nwc,mesh_loc%me))
1643  msup = maxval(m_loc_to_glob)
1644  minf = minval(m_loc_to_glob)
1645  DO m = 1, mesh_loc%me
1646  DO n = 1, nwc
1647  mop = mesh%neigh(n,m_loc_to_glob(m))
1648  IF (mop==0) THEN
1649  mesh_loc%neigh(n,m) = 0
1650  ELSE IF(mop<minf .OR. mop>msup) THEN
1651  !mesh_loc%neigh(n,m) = mop
1652  mesh_loc%neigh(n,m) = -mop !JLG AUG 13 2015
1653  ELSE
1654  mesh_loc%neigh(n,m) = m_glob_to_loc(mop)
1655  END IF
1656  END DO
1657  END DO
1658  !==End re-order neigh
1659 
1660  !==Re-order i_d
1661  ALLOCATE(mesh_loc%i_d(mesh_loc%me))
1662  mesh_loc%i_d = mesh%i_d(m_loc_to_glob)
1663  !==End re-order i_d
1664 
1665  !==Re-order neighs
1666  ALLOCATE(mesh_loc%neighs(mesh_loc%mes))
1667  mesh_loc%neighs = m_glob_to_loc(mesh%neighs(mes_loc(1):mes_loc(2)))
1668  !==End re-order neighs
1669 
1670 
1671  !==Re-order sides
1672  ALLOCATE(mesh_loc%sides(mesh_loc%mes))
1673  mesh_loc%sides = mesh%sides(mes_loc(1):mes_loc(2))
1674  !==End re-order sides
1675 
1676  !==Re-order jjs
1677  ALLOCATE(mesh_loc%jjs(nws,mesh_loc%mes))
1678  DO ns = 1, nws
1679  mesh_loc%jjs(ns,:) = glob_to_loc(mesh%jjs(ns,mes_loc(1):mes_loc(2)))
1680  END DO
1681  !==End re-order jjs
1682 
1683  !TEST
1684 
1685  DO ms = 1, mesh_loc%mes
1686  m = mesh_loc%neighs(ms)
1687  DO ns = 1, nws
1688  test = .true.
1689  DO n = 1, nw
1690  IF (maxval(abs(mesh_loc%rr(:,mesh_loc%jj(n,m))-mesh_loc%rr(:,mesh_loc%jjs(ns,ms)))) .LT. 1.d-10) THEN
1691  test = .false.
1692  END IF
1693  END DO
1694  IF (test) THEN
1695  WRITE(*,*) 'bug in create local mesh, non consistent numbering'
1696  END IF
1697  END DO
1698  END DO
1699 
1700 
1701  !TEST
1702 
1703  DEALLOCATE(m_loc_to_glob)
1704 
1705  END SUBROUTINE create_local_mesh
1706 
1707  SUBROUTINE free_mesh(mesh)
1708  USE def_type_mesh
1709  USE my_util
1710  IMPLICIT NONE
1711  TYPE(mesh_type) :: mesh
1712 
1713  DEALLOCATE(mesh%jj)
1714  DEALLOCATE(mesh%jjs)
1715  DEALLOCATE(mesh%rr)
1716  DEALLOCATE(mesh%neigh)
1717  DEALLOCATE(mesh%sides)
1718  DEALLOCATE(mesh%neighs)
1719  DEALLOCATE(mesh%i_d)
1720 
1721  NULLIFY(mesh%loc_to_glob)
1722  NULLIFY(mesh%disp)
1723  NULLIFY(mesh%domnp)
1724  NULLIFY(mesh%j_s)
1725 
1726  IF (mesh%edge_stab) THEN
1727  DEALLOCATE(mesh%iis)
1728  NULLIFY(mesh%jji)
1729  DEALLOCATE(mesh%jjsi)
1730  DEALLOCATE(mesh%neighi)
1731  END IF
1732 
1733  mesh%dom_me = 0
1734  mesh%dom_np = 0
1735  mesh%dom_mes = 0
1736  mesh%me = 0
1737  mesh%mes = 0
1738  mesh%np = 0
1739  mesh%nps = 0
1740  mesh%mi =0
1741  mesh%edge_stab = .false.
1742 
1743  END SUBROUTINE free_mesh
1744 
1745  SUBROUTINE free_interface(interf)
1746  USE def_type_mesh
1747  USE my_util
1748  IMPLICIT NONE
1749  TYPE(interface_type) :: interf
1750 
1751  interf%mes = 0
1752  DEALLOCATE(interf%mesh1)
1753  DEALLOCATE(interf%mesh2)
1754  DEALLOCATE(interf%jjs1)
1755  DEALLOCATE(interf%jjs2)
1756  END SUBROUTINE free_interface
1757 
1758  SUBROUTINE reassign_per_pts(mesh, partition, list_pts)
1759  USE def_type_mesh
1760  USE my_util
1761  IMPLICIT NONE
1762 
1763  TYPE(mesh_type) , INTENT(IN) :: mesh
1764  INTEGER, DIMENSION(mesh%me), INTENT(INOUT) :: partition
1765  INTEGER, DIMENSION(mesh%np,3), INTENT(IN) :: list_pts
1766 
1767  INTEGER :: i, j_loc, proc_min, index, i_loc, m, mop, n, proc1, proc2
1768  INTEGER, DIMENSION(50) :: list_elmts
1769  LOGICAL :: okay
1770 
1771 
1772  list_elmts = 0
1773  DO i = 1, mesh%np
1774  IF (list_pts(i,2)==0) cycle
1775  j_loc = list_pts(i,1)
1776  list_elmts=0
1777  index = 1
1778  list_elmts(index) = list_pts(i,2)
1779  okay = .true.
1780  DO WHILE (okay)
1781  m=list_elmts(index)
1782  okay = .false.
1783  i_loc=index
1784  DO n = 1, 3
1785  mop = mesh%neigh(n, m)
1786  IF (mop == 0) cycle
1787  IF (minval(abs(mesh%jj(:,mop)-j_loc)) /=0) cycle
1788  IF (minval(abs(mop-list_elmts))==0) cycle
1789  okay=.true.
1790  i_loc = i_loc + 1
1791  IF (i_loc-index==2) THEN
1792  CALL error_petsc('BUG in reassign_per_pts, how is that possible?')
1793  END IF
1794  list_elmts(i_loc) = mop
1795  END DO
1796  index = i_loc
1797  END DO
1798 !!$ WRITE(*,*) i, list_elmts(1:index)
1799  IF (list_pts(i,3) == 0) THEN ! point au bord du bord periodique, ou sur une arete
1800  proc_min = minval(partition(list_elmts(1:index)))
1801  partition(list_elmts(1)) = proc_min
1802  ELSE ! deux elements du bord periodique touchent le point
1803  IF (list_elmts(index) /= list_pts(i,3)) THEN
1804  CALL error_petsc('BUG in reassign_per_pts, wrong element')
1805  END IF
1806  proc1 = partition(list_elmts(1))
1807  proc2 = partition(list_elmts(2))
1808  partition(list_elmts(2:index-1)) = min(proc1,proc2)
1809  END IF
1810  END DO
1811 
1812  END SUBROUTINE reassign_per_pts
1813 
1814 
1815 END MODULE metis_sfemans
def_type_mesh::interface_type
Definition: def_type_mesh.f90:116
sub_plot
Definition: sub_plot.f90:4
def_type_mesh::periodic_data
Definition: def_type_mesh.f90:12
metis_sfemans
Definition: metis_distribution_sfemans.f90:1
def_type_mesh::mesh_type
Definition: def_type_mesh.f90:74
my_util
Definition: my_util.f90:1
my_util::error_petsc
subroutine error_petsc(string)
Definition: my_util.f90:16
metis_sfemans::free_interface
subroutine, public free_interface(interf)
Definition: metis_distribution_sfemans.f90:1746
sub_plot::plot_const_p1_label
subroutine plot_const_p1_label(jj, rr, uu, file_name)
Definition: sub_plot.f90:265
metis_sfemans::reassign_per_pts
subroutine reassign_per_pts(mesh, partition, list_pts)
Definition: metis_distribution_sfemans.f90:1759
metis_sfemans::part_mesh_m_t_h_phi
subroutine, public part_mesh_m_t_h_phi(nb_proc, list_u, list_T_in, list_h_in, list_phi, mesh, list_of_interfaces, part, my_periodic)
Definition: metis_distribution_sfemans.f90:14
metis_sfemans::free_mesh
subroutine, public free_mesh(mesh)
Definition: metis_distribution_sfemans.f90:1708
metis_sfemans::metis_noptions
integer metis_noptions
Definition: metis_distribution_sfemans.f90:8
metis_sfemans::part_mesh_mhd
subroutine, public part_mesh_mhd(nb_proc, vwgt, mesh, list_of_interfaces, part, my_periodic)
Definition: metis_distribution_sfemans.f90:552
metis_sfemans::metis_option_numbering
integer metis_option_numbering
Definition: metis_distribution_sfemans.f90:8
metis_sfemans::part_mesh_mhd_bis
subroutine, public part_mesh_mhd_bis(nb_proc, list_u, list_h_in, list_phi, mesh, list_of_interfaces, part, my_periodic)
Definition: metis_distribution_sfemans.f90:834
metis_sfemans::extract_mesh
subroutine, public extract_mesh(communicator, nb_proc, mesh_glob, part, list_dom, mesh, mesh_loc)
Definition: metis_distribution_sfemans.f90:1327
metis_sfemans::create_local_mesh
subroutine create_local_mesh(mesh, mesh_loc, me_loc, mes_loc, np_loc)
Definition: metis_distribution_sfemans.f90:1505
def_type_mesh
Definition: def_type_mesh.f90:4
metis_sfemans::epsilon
real(kind=8) epsilon
Definition: metis_distribution_sfemans.f90:6