GetFEM  5.4.2
getfem_enumeration_dof_para.cc
1 /*===========================================================================
2 
3  Copyright (C) 2012-2020 Yves Renard, Julien Pommier.
4 
5  This file is a part of GetFEM
6 
7  GetFEM is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program; if not, write to the Free Software Foundation,
18  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
19 
20 ===========================================================================*/
21 #include <queue>
22 #include <queue>
23 #include "getfem/dal_singleton.h"
24 #include "getfem/getfem_mesh_fem.h"
25 
26 
27 namespace getfem {
28 
29 #if GETFEM_PARA_LEVEL > 1
30 
31  struct fem_dof {
32  size_type ind_node;
33  pdof_description pnd;
34  size_type part;
35  };
36 
37 /* Fonction de renumérotation des degrés de libertés */
38 
39 /// Parallel Enumeration of dofs
40 void mesh_fem::enumerate_dof_para(void) const {
41 
42 #if 0
43 
44  GMM_TRACE2("Enumeration dof para !!!!!!!!!!!!!!!!");
45  GMM_ASSERT1(linked_mesh_ != 0, "Uninitialized mesh_fem");
46  context_check();
47  if (fe_convex.card() == 0) {
48  dof_enumeration_made = true;
49  nb_total_dof = 0;
50  return;
51  }
52 
53  dof_structure.clear();
54  MPI_Status mstatus;
55  fem_dof fd;
56 
57 /* Récupération du nombre total de régions (procs)!!!!*/
58  int num_rg, nb_rg;
59  MPI_Comm_rank(MPI_COMM_WORLD, &num_rg);
60  MPI_Comm_size(MPI_COMM_WORLD, &nb_rg);
61 
62 
63  cout<<"nb total region : "<<nb_rg<<" et nb_points = "<<linked_mesh().nb_points()<<endl;
64 
65 /* Récupération du num de la région */
66  //size_type num_rg;
67  //num_rg = linked_mesh.mpi_region.id();
68 
69 /* Création de la liste des ddl */
70 /// list_of_dof[1:nb_total_dof_mesh, 1:2]
71  std::vector<size_type> list_of_dof_numeration;
72  std::vector<size_type> list_of_dof_num_rg;
73  dal::bit_vector enumeration_of_dof_made;
74 
75  list_of_dof_numeration.resize(100000);
76  list_of_dof_num_rg.resize(100000);
77  // enumeration_of_dof_made(100000);
78 
79 /* Création de la liste des ddl interface */
80  std::vector<size_type> list_of_dof_linkable_index;
81  std::vector<size_type> list_of_dof_linkable_to;
82 
83  list_of_dof_linkable_index.resize(10000);
84  list_of_dof_linkable_to.resize(10000);
85 
86 /* Création de la liste des ddl globaux */
87  std::vector<size_type> list_of_global_dof_index;
88  std::vector<size_type> list_of_global_dof_local_index;
89 
90  list_of_global_dof_index.resize(10000);
91  list_of_global_dof_local_index.resize(10000);
92 
93 /* Initialisation de l'itérateur sur list_of_dof */
94  size_type iter_dof = 0;
95 
96 /* Construction de la liste des cv à la charge de chaque proc en fonction de la région */
97 /// list_of_cv[1:nb_cv_total_mesh, 1:2]
98  std::vector<size_type> list_of_cv_num_rg;
99  //list_of_cv_num_rg = size_type(0);
100  std::vector<size_type> list_of_cv_first_index;
101 
102  std::vector<size_type> list_of_icv;
103  std::vector<size_type> icv_in_list;
104 
105  list_of_cv_num_rg.resize(10000);
106  list_of_cv_first_index.resize(10000);
107  list_of_icv.resize(1000);
108  icv_in_list.resize(1000);
109 
110 /// En parallèle sur les régions on récupère les indices des cv de chaque proc
111  //dal::bit_vector index_tab;
112  const std::vector<size_type> &cmk = linked_mesh().cuthill_mckee_ordering();
113  //index_tab = linked_mesh().region(num_rg).index();
114 
115  cout<<"cmk.size = "<<cmk.size()<<endl;
116  cout<<"cmk[0] = "<<cmk[0]<<endl;
117  cout<<"cml[10] = "<<cmk[10]<<endl;
118  cout<<"nb_cv = "<<linked_mesh().convex_index().card()<<endl;
119 
120  size_type nb_cv = 0;
121  std::vector<size_type> neighbors;
122  bgeot::pgeotrans_precomp pgp = 0;
123  base_node P;
124  bgeot::pgeotrans_precomp pgpj = 0;
125  base_node Pj;
126 
127 // Boucle i pour remplir la liste des cv:
128  GMM_TRACE2("Initialisation of lists cv");
129  // for(dal::bv_visitor i(index_tab); !i.finished(); ++i)
130  cout<<"me = "<<num_rg<<endl;
131 
132  bool entre = false;
133 
134  cout<<"bool avant : "<<entre<<endl;
135 
136  for(size_type i = cmk[0]; i<cmk.size(); i++)
137  {
138  size_type icv = cmk[i];
139  if(linked_mesh().region(num_rg).is_in(icv))
140  {
141  GMM_TRACE2("ICI 0");
142  // cout<<"i = "<<i<<endl;
143  //cout<<"me "<<num_rg<<" et icv = "<<icv<<endl;
144  list_of_cv_num_rg[nb_cv] = num_rg;
145  GMM_TRACE2("ICI 1");
146  list_of_cv_first_index[nb_cv] = iter_dof;
147  GMM_TRACE2("ICI 2");
148  list_of_icv[nb_cv] = icv;
149  icv_in_list[icv] = nb_cv;
150 
151  pfem pf = fem_of_element(icv);
152  size_type nbd = pf->nb_dof(icv);
153  iter_dof += nbd;
154  nb_cv += 1;
155  //cout<<"la!!!!!!"<<nb_cv<<endl;
156  entre = true;
157  }
158  else
159  {
160  list_of_cv_num_rg[nb_cv] = 0;
161  list_of_cv_first_index[nb_cv] = 0;
162  list_of_icv[nb_cv]=0;
163  icv_in_list[icv]=0;
164 
165  pfem pf = fem_of_element(icv);
166  size_type nbd = pf->nb_dof(icv);
167  iter_dof += nbd;
168  nb_cv += 1;
169  //cout<<"ou la !!!!!!!!!"<<nb_cv<<endl;
170 
171  }
172  //cout<<"me = "<<num_rg<<"iteration i : "<<i<<endl;
173  }
174 
175  cout<<"bool : "<<entre<<endl;
176 
177  //int nb_cv_tot;
178 // Mise en commun par échange MPI_AllReduce du nombre totale de cv sur le mesh
179  //GMM_TRACE2("Echange MPI 1");
180  // MPI_Allreduce(&nb_cv, &nb_cv_tot, 1, MPI_INTEGER, MPI_SUM, MPI_COMM_WORLD);
181 
182  cout<<"nb_cv_tot = "<<nb_cv<<endl;
183 
184 
185 
186 
187  std::vector<size_type> list_of_cv_num_rg_Recv;
188  std::vector<size_type> list_of_cv_first_index_Recv;
189  std::vector<size_type> list_of_icv_Recv;
190  std::vector<size_type> icv_in_list_Recv;
191 
192  cout<<"size(cv_num_rg) = "<<list_of_cv_num_rg.size()<<endl;
193 
194  list_of_cv_num_rg.resize(nb_cv);
195  list_of_cv_num_rg_Recv.resize(nb_cv);
196  list_of_cv_first_index_Recv.resize(nb_cv);
197  list_of_icv_Recv.resize(nb_cv);
198  icv_in_list_Recv.resize(nb_cv);
199 
200  MPI_Barrier(MPI_COMM_WORLD);
201 
202 // Mise en commun par échange MPI_AllReduce de la liste list_of_cv_num_rg
203  GMM_TRACE2("Echange MPI 2");
204  MPI_Allreduce (&list_of_cv_num_rg[0], &list_of_cv_num_rg_Recv, nb_cv,
205  MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD);
206 
207  cout<<"num_rg[10] = "<<list_of_cv_num_rg[10]<<endl;
208  cout<<"num_rg_Recv[10] = "<<list_of_cv_num_rg_Recv[10]<<endl;
209 
210  GMM_TRACE2("Echange 3");
211 // Mise en commun par échange MPI_AllReduce de la liste list_of_cv_num_rg
212  MPI_Allreduce (&list_of_cv_first_index[0], &list_of_cv_first_index_Recv, nb_cv,
213  MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD);
214 
215  GMM_TRACE2("Echange 4");
216  MPI_Allreduce (&list_of_icv[0], &list_of_icv_Recv, nb_cv,
217  MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD);
218 
219  GMM_TRACE2("Echange 5");
220 
221  MPI_Allreduce (&icv_in_list[0], &icv_in_list_Recv, nb_cv,
222  MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD);
223 
224 /* Construction de la liste des ddl à la charge de chaque proc */
225  size_type nb_dof_rg = 0;
226  size_type nb_dof_tot;
227  //size_type nb_dof_inter = 0;
228  size_type nb_global_dof = 0;
229  size_type nb_global_dof_tot;
230  size_type nb_dof_linkable = 0;
231 
232 // Pour chaque cv dont ce proc a la charge :
233  GMM_TRACE2("Attributions des ddl aux rg");
234  for(size_type cv = 0; cv < list_of_cv_num_rg_Recv.size(); ++cv)
235  {
236  size_type icv = list_of_icv_Recv[cv];
237  if (list_of_cv_num_rg_Recv[cv] == num_rg)
238  {
239  pfem pf = fem_of_element(icv);
240  size_type nbd = pf->nb_dof(icv);
241  nb_dof_rg += nbd;
242  pdof_description andof = global_dof(pf->dim());
243 
244 // pour chaque ddl associé à ce cv :
245  for (size_type i = list_of_cv_first_index_Recv[cv];
246  i <= list_of_cv_first_index_Recv[cv] + nbd; i++)
247  {
248  fd.pnd = pf->dof_types()[i];
249  fd.part = get_dof_partition(icv);
250 
251 // Test si le ddl est raccordable
252  if (dof_linkable(fd.pnd))
253  {
254  size_type bool_rg = 0;
255  size_type bool_inter = 0;
256  P.resize(linked_mesh().dim());
257  pgp->transform(linked_mesh().points_of_convex(icv), i, P);
258 
259 // Récupération des voisins qui possèdent ce même ddl :
260  neighbors = linked_mesh().convex_to_point(i);
261  for (size_type jcv = neighbors[0]; jcv < neighbors.size(); ++jcv)
262  {
263 // Si le voisin appartient à la même région (ie pas ddl interface)
264  if (list_of_cv_num_rg_Recv[icv_in_list_Recv[jcv]] == num_rg)
265  {
266  bool_rg++;
267  }
268 // Sinon si c'est un dof interface "et" qui doit être à la charge de cette région
269  else if (list_of_cv_num_rg_Recv[icv_in_list_Recv[jcv]] > num_rg)
270  {
271  bool_inter++;
272  }
273  }
274 // Test si pas ddl interface
275  if (bool_rg == neighbors.size() || bool_inter+bool_rg == neighbors.size())
276  // ie tout les voisins raccordable sont dans cette même region
277  {
278  /* for(size_type jcv = neighbors[0]; jcv < neighbors.size(); ++jcv)
279  {
280  list_of_dof_linkable_index[nb_dof_linkable] = list_of_cv_first_index_Recv[jcv]+j;
281  list_of_dof_linkable_to[nb_dof_linkable] = i;
282  nb_dof_linkable ++;
283  }
284  list_of_dof_num_rg[i] = num_rg;
285  }
286  else if ((bool_inter + bool_rg)==neighbors.size())
287  // ie tout les voisins raccordable doivent être associé à cette region
288  {*/
289  list_of_dof_num_rg[i] = num_rg;
290  for (size_type jcv = neighbors[0]; jcv < neighbors.size(); ++jcv)
291  {
292 /// on associe le ddl correspondant au proc
293  pfem pfj = fem_of_element(jcv);
294  size_type nbdj = pfj->nb_dof(jcv);
295  for(size_type j = 0; j < nbdj; j++)
296  {
297  Pj.resize(linked_mesh().dim());
298  pgpj->transform(linked_mesh().points_of_convex(jcv), j, Pj);
299  if (&P == &Pj)
300  {
301  list_of_dof_linkable_index[nb_dof_linkable] =
302  list_of_cv_first_index_Recv[icv_in_list_Recv[jcv]]+j;
303  list_of_dof_linkable_to[nb_dof_linkable] = i;
304  nb_dof_linkable++;
305  list_of_dof_num_rg[list_of_cv_first_index_Recv
306  [icv_in_list_Recv[jcv]]+j] = num_rg;
307  }
308  }
309  }
310  }
311 
312  }
313 // Si ddl global
314  else if(fd.pnd == andof)
315  {
316 /// on recupère son numéro global : encountered_global_dof[1:3, num_region] = [num_global, icv, i]
317  size_type num = pf->index_of_global_dof(icv, i);
318  list_of_global_dof_index[nb_global_dof] = num;
319  list_of_global_dof_local_index [nb_global_dof] = i;
320  list_of_dof_num_rg[i] = num_rg;
321  nb_global_dof++;
322  }
323 // Si le ddl est non raccordable
324  else
325  {
326 /// on associe ce ddl au proc => list_of_dof[i, 1] = num_region (ou qqch de remarquable!!!!)
327  list_of_dof_num_rg[i] = num_rg;
328  } //end if
329  } // end for
330  } // end if
331  } // end for
332 
333 
334 
335 // Mise en commun par échange MPI_AllReduce de la liste list_of _dof
336 
337 // Mise en commun par échange MPI_AllReduce du nombre totale de dof sur le mesh
338 
339  MPI_Allreduce (&nb_dof_rg, &nb_dof_tot, 1, MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD);
340  size_type nbd_p;
341  size_type numerot = 0;
342  for (int p = 0; p < nb_rg; p++)
343  {
344  // Si il a des proc plus petit que num_rg, il recoit le nb de dof des autres plus petit pour mettre à jour son indice de début de numérotation
345  if (p < num_rg)
346  {
347  MPI_Recv(&nbd_p, 1, MPI_UNSIGNED, p, 100, MPI_COMM_WORLD, &mstatus);
348  numerot += nbd_p;
349  }
350  // Sinon il envoi le nombre de dof qu'il a à sa charge au autre qui lui sont supérieur
351  else if (p > num_rg)
352  {
353  MPI_Send(&nb_dof_rg, 1, MPI_UNSIGNED, p, 100, MPI_COMM_WORLD);
354  }
355  }
356 
357 
358  std::vector<size_type> list_of_dof_num_rg_Recv;
359  list_of_dof_num_rg_Recv.resize(nb_dof_tot);
360 
361 // Mise en commun par échange MPI_AllReduce de la liste list_of_cv_num_rg
362  MPI_Allreduce(&list_of_dof_num_rg[0], &list_of_dof_num_rg_Recv, nb_dof_tot,
363  MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD);
364 
365 
366  std::vector<size_type> list_of_global_dof_index_Recv;
367  std::vector<size_type> list_of_global_dof_local_index_Recv;
368  size_type nb_global_dof_Recv;
369 
370 // Mise en commun des information sur les global_dof
371  MPI_Allreduce (&nb_global_dof, &nb_global_dof_tot, 1, MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD);
372 
373  list_of_global_dof_index_Recv.resize(nb_global_dof_tot);
374  list_of_global_dof_local_index_Recv.resize(nb_global_dof_tot);
375 
376  for (int p = 0; p<nb_rg; p++)
377  {
378  MPI_Send (&nb_global_dof, 1, MPI_UNSIGNED, p, 200, MPI_COMM_WORLD);
379  MPI_Recv (&nb_global_dof_Recv, 1, MPI_UNSIGNED, p, 200, MPI_COMM_WORLD, &mstatus);
380 
381  MPI_Send(&list_of_global_dof_index[0], nb_global_dof, MPI_UNSIGNED, p, 300, MPI_COMM_WORLD);
382 
383  MPI_Recv (&list_of_global_dof_index_Recv[0], nb_global_dof_Recv,
384  MPI_UNSIGNED, p, 300, MPI_COMM_WORLD, &mstatus);
385 
386  MPI_Send(&list_of_global_dof_local_index[0], nb_global_dof, MPI_UNSIGNED, p, 400, MPI_COMM_WORLD);
387 
388  MPI_Recv (&list_of_global_dof_local_index_Recv[0], nb_global_dof_Recv,
389  MPI_UNSIGNED, p, 400, MPI_COMM_WORLD, &mstatus);
390  }
391 
392 
393 
394 /* Numérotation des ddl en charge de chaque processeur */
395 // Pour chaque cv dont ce proc a la charge :
396  GMM_TRACE2("Numerotation des ddl");
397  size_type ind_linkable = 0;
398  for(size_type icv = 0; icv < list_of_cv_num_rg_Recv.size(); ++icv)
399  {
400  pfem pf = fem_of_element(icv);
401  size_type nbd = pf->nb_dof(icv);
402 
403 // pour chaque ddl associé à ce cv :
404  for (size_type i = list_of_cv_first_index_Recv[icv];
405  i <= list_of_cv_first_index_Recv[icv] + nbd; i++)
406  {
407  if (list_of_dof_num_rg_Recv[i] == num_rg && !enumeration_of_dof_made[i])
408  {
409  if (!dof_linkable(fd.pnd))
410  {
411  list_of_dof_numeration[i] = numerot;
412  numerot += Qdim / pf->target_dim();
413  }
414 // Test si ddl raccordable
415  else if (dof_linkable(fd.pnd))
416  {
417 /// Recherche des cv ayant ce même point et dont le proc à la charge
418  while (list_of_dof_linkable_to[ind_linkable] == i)
419  {
420 /// Récupération des indices correspondants dans list_of_dof et Numérotation dans list_of_dof(indices)
421  list_of_dof_numeration[list_of_dof_linkable_index[ind_linkable]] = numerot;
422  enumeration_of_dof_made[list_of_dof_linkable_index[ind_linkable]] = true;
423  ind_linkable++;
424  }
425  list_of_dof_numeration[i] = numerot;
426  enumeration_of_dof_made[i] = true;
427  numerot += Qdim / pf->target_dim();
428  } // Fin Si
429  } // Fin boucle
430  } // Fin Boucle
431  }// Fin Boucle
432 
433 // Traitement des ddl globaux
434 // Boucle sur list_of_global_dof_in_charge
435 /* if(num_rg == 0)// temporairement
436  {
437  for (size_type i=0; i < list_of_global_dof_index_Recv.size(); i++)
438  {
439  pfem pf = fem_of_element(icv);
440 
441  if(!enumeration_of_dof_made[i])
442  {
443 /// Récupère les indices ayant le même num_global
444  for (size_type j = i; j < list_of_global_dof_index_Recv.size(); j++)
445  {
446 /// Numérotation
447  if (list_of_global_dof_index_Recv[j] == list_of_global_dof_index_Recv[i]
448  && !enumeration_of_dof_made[j])
449  {
450  list_of_dof_numeration[list_of_global_dof_local_index_Recv[j]] = numerot;
451  enumeration_of_dof_made[list_of_gloabl_dof_local_index_Recv[j]] = true;
452  }
453  }
454  list_of_dof_numeration[list_of_global_dof_local_index[i]] = numerot;
455  enumeration_of_dof_made[list_of_global_dof_local_index[i]] = true;
456  numerot += Qdim / pf->target_dim();
457  }
458  }
459  }*/
460 // Fin boucle
461 
462 // Mise en commun de list_of_dof par échange avec MPI_AllReduce
463  std::vector<size_type> list_of_dof_numeration_Recv;
464  list_of_dof_numeration_Recv.resize(nb_dof_tot);
465 
466  MPI_Allreduce(&list_of_dof_numeration[0], &list_of_dof_numeration_Recv, numerot,
467  MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD);
468 
469 // Envoi de la structure numérotée
470  GMM_TRACE2("Envoi de la numerotation");
471  std::vector<size_type> tab;
472  size_type ind_tab = 0;
473  for(size_type icv = 0; icv < list_of_cv_num_rg.size(); ++icv)
474  {
475  pfem pf = fem_of_element(icv);
476  if (list_of_cv_num_rg[icv] == num_rg)
477  {
478  size_type nbd = pf->nb_dof(icv);
479  tab.resize(nbd);
480  for (size_type i = list_of_cv_first_index_Recv[icv]; i < list_of_cv_first_index_Recv[icv] + nbd; i++)
481  {
482  tab[ind_tab] = list_of_dof_numeration[i];
483  ind_tab++;
484  }
485  }
486  dof_structure.add_convex_noverif(pf->structure(icv), tab.begin(), icv);
487  }
488 
489  nb_total_dof = nb_dof_tot;
490 
491 
492 #endif
493 }
494 
495 #endif
496 
497 
498 } // end of getfem namespace
getfem::pdof_description
dof_description * pdof_description
Type representing a pointer on a dof_description.
Definition: getfem_fem.h:160
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
dal_singleton.h
A simple singleton implementation.
getfem::mesh_fem::fem_of_element
virtual pfem fem_of_element(size_type cv) const
Return the basic fem associated with an element (if no fem is associated, the function will crash!...
Definition: getfem_mesh_fem.h:444
bgeot::mesh_structure::clear
void clear()
erase the mesh
Definition: bgeot_mesh_structure.cc:205
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
bgeot::mesh_structure::convex_index
const dal::bit_vector & convex_index() const
Return the list of valid convex IDs.
Definition: bgeot_mesh_structure.h:91
getfem_mesh_fem.h
Define the getfem::mesh_fem class.
getfem::dof_linkable
bool dof_linkable(pdof_description)
Says if the dof is linkable.
Definition: getfem_fem.cc:615
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:244
getfem::mesh_fem::linked_mesh
const mesh & linked_mesh() const
Return a reference to the underlying mesh.
Definition: getfem_mesh_fem.h:279
bgeot::mesh_structure::convex_to_point
const ind_cv_ct & convex_to_point(size_type ip) const
Return a container of the convexes attached to point ip.
Definition: bgeot_mesh_structure.h:166
getfem::mesh::nb_points
size_type nb_points() const
Give the number of geometrical nodes in the mesh.
Definition: getfem_mesh.h:194
getfem::context_dependencies::context_check
bool context_check() const
return true if update_from_context was called
Definition: getfem_context.h:126
getfem::mesh::cuthill_mckee_ordering
const std::vector< size_type > & cuthill_mckee_ordering() const
Return the list of convex IDs for a Cuthill-McKee ordering.
getfem::global_dof
pdof_description global_dof(dim_type d)
Description of a global dof, i.e.
Definition: getfem_fem.cc:542