GetFEM++  5.3
bgeot_mesh_structure.cc
1 /*===========================================================================
2 
3  Copyright (C) 1999-2017 Yves Renard
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 
22 
23 
25 
26 namespace bgeot {
27 
28 
29  dal::bit_vector mesh_structure::convex_index(dim_type n) const {
30  dal::bit_vector res = convex_tab.index();
31  for (dal::bv_visitor cv(convex_tab.index()); !cv.finished(); ++cv)
32  if (structure_of_convex(cv)->dim() != n) res.sup(cv);
33  return res;
34  }
35 
37  size_type ip) const {
38  mesh_convex_structure::ind_pt_ct::const_iterator it;
39  size_type ind = 0;
40  for (it=convex_tab[ic].pts.begin();
41  it != convex_tab[ic].pts.end() && (*it) != ip; ++it) ++ind;
42  GMM_ASSERT1(it != convex_tab[ic].pts.end(),
43  "This point does not exist on this convex.");
44  return ind;
45  }
46 
48  if (i == j) return;
49  std::vector<size_type> doubles;
50 
51  for (size_type k = 0; k < points_tab[i].size(); ++k) {
52  size_type cv = points_tab[i][k];
53  for (size_type l = 0; l < convex_tab[cv].pts.size(); ++l) {
54  size_type &ind = convex_tab[cv].pts[l];
55  if (ind == i) ind = j;
56  else if (ind == j) { ind = i; doubles.push_back(cv); }
57  }
58  }
59  for (size_type k = 0; k < points_tab[j].size(); ++k) {
60  size_type cv = points_tab[j][k];
61  if (std::find(doubles.begin(), doubles.end(), cv) == doubles.end()) {
62  for (size_type l = 0; l < convex_tab[cv].pts.size(); ++l)
63  if (convex_tab[cv].pts[l] == j) convex_tab[cv].pts[l] = i;
64  }
65  }
66  points_tab.swap(i,j);
67  }
68 
70  if (i == j) return;
71  std::vector<size_type> doubles;
72 
73  if (is_convex_valid(i))
74  for (size_type k = 0; k < convex_tab[i].pts.size(); ++k) {
75  size_type ip = convex_tab[i].pts[k];
76  for (size_type l = 0; l < points_tab[ip].size(); ++l) {
77  size_type &ind = points_tab[ip][l];
78  if (ind == i) ind = j;
79  else if (ind == j) { ind = i; doubles.push_back(ip); }
80  }
81  }
82  if (is_convex_valid(j))
83  for (size_type k = 0; k < convex_tab[j].pts.size(); ++k) {
84  size_type ip = convex_tab[j].pts[k];
85  if (std::find(doubles.begin(), doubles.end(), ip) == doubles.end()) {
86  for (size_type l = 0; l < points_tab[ip].size(); ++l)
87  if (points_tab[ip][l] == j) points_tab[ip][l] = i;
88  }
89  }
90  convex_tab.swap(i,j);
91  }
92 
93  size_type mesh_structure::add_segment(size_type a, size_type b) {
94  static pconvex_structure cs = NULL;
95  if (!cs) cs = simplex_structure(1);
96  size_type t[2]; t[0] = a; t[1] = b;
97  return add_convex(cs, &t[0]);
98  }
99 
101  if (!(is_convex_valid(ic))) return;
102  for (size_type l = 0; l < convex_tab[ic].pts.size(); ++l) {
103  size_type &ind = convex_tab[ic].pts[l];
104  std::vector<size_type>::iterator it1= points_tab[ind].begin(), it2 = it1;
105  std::vector<size_type>::iterator ite= points_tab[ind].end();
106  for (; it2 != ite; ++it2) { *it1 = *it2; if (*it1 != ic) ++it1; }
107  points_tab[ind].pop_back();
108  }
109  convex_tab.sup(ic);
110  }
111 
113  return add_convex( (structure_of_convex(ic)->faces_structure())[f],
114  ind_points_of_face_of_convex(ic, f).begin());
115  }
116 
119  for (short_type iff = 0; iff < ps->nb_faces(); ++iff)
120  add_convex( (ps->faces_structure())[iff],
121  ind_points_of_face_of_convex(ic, iff).begin());
122  }
123 
124  void mesh_structure::to_faces(dim_type n) { // not efficient
125  dal::bit_vector nn = convex_index();
126  for (dal::bv_visitor i(nn); !i.finished(); ++i)
127  if (convex_tab[i].cstruct->dim() == n)
128  { add_faces_of_convex(i); sup_convex(i); }
129  }
130 
131  void mesh_structure::to_edges(void) { // not efficient at all !
132  dim_type dmax = 0;
133  dal::bit_vector nn = convex_index();
134  for (dal::bv_visitor i(nn); !i.finished(); ++i)
135  dmax = std::max(dmax, convex_tab[i].cstruct->dim());
136  for ( ; dmax > 1; --dmax) to_faces(dmax);
137  }
138 
139  size_type mesh_structure::nb_convex_with_edge(size_type i1, size_type i2) {
140  size_type nb = 0;
141  for (size_type k = 0; k < points_tab[i1].size(); ++k) {
142  size_type cv = points_tab[i1][k];
143  for (size_type l = 0; l < convex_tab[cv].pts.size(); ++l)
144  if (convex_tab[cv].pts[l] == i2) { ++nb; break; }
145  }
146  return nb;
147  }
148 
149  void mesh_structure::convex_with_edge(size_type i1, size_type i2,
150  std::vector<size_type> &ipt) const {
151  ipt.resize(0);
152  for (size_type k = 0; k < points_tab[i1].size(); ++k) {
153  size_type cv = points_tab[i1][k];
154  for (size_type l = 0; l < convex_tab[cv].pts.size(); ++l)
155  if (convex_tab[cv].pts[l] == i2) { ipt.push_back(cv); break; }
156  }
157  }
158 
159  void mesh_structure::ind_points_to_point(size_type ip, ind_set &s) const {
160  // not efficient
161  s.resize(0);
162  for (size_type k = 0; k < points_tab[ip].size(); ++k) {
163  size_type cv = points_tab[ip][k];
164  for (size_type l = 0; l < convex_tab[cv].pts.size(); ++l)
165  if (ip != convex_tab[cv].pts[l]) {
166  size_type ind = convex_tab[cv].pts[l];
167  if (std::find(s.begin(), s.end(), ind) != s.end()) s.push_back(ind);
168  }
169  }
170  }
171 
174  short_type iff) const {
175  const mesh_convex_structure *q = &(convex_tab[ic]);
176  std::vector<size_type>::const_iterator r = q->pts.begin();
177  const convex_ind_ct *p = &(q->cstruct->ind_points_of_face(iff));
178  return ind_pt_face_ct(r, p->begin(), p->end());
179  }
180 
181  size_type mesh_structure::memsize(void) const {
182  size_type mems = sizeof(mesh_structure) + points_tab.memsize()
183  + convex_tab.memsize();
184  for (size_type i = 0; i < convex_tab.size(); ++i)
185  mems += convex_tab[i].pts.size() * sizeof(size_type);
186  for (size_type i = 0; i < points_tab.size(); ++i)
187  mems += points_tab[i].size() * sizeof(size_type);
188  return mems;
189  }
190 
192  size_type i, j = nb_convex();
193  for (i = 0; i < j; i++)
194  if (!convex_tab.index_valid(i))
195  swap_convex(i, convex_tab.ind_last());
196 
197  if (points_tab.size())
198  for (i=0, j = (points_tab.end()-points_tab.begin())-1; i < j; ++i, --j){
199  while (i < j && !(points_tab[i].empty())) ++i;
200  while (i < j && points_tab[j].empty()) --j;
201  if (i < j) swap_points(i, j);
202  }
203  }
204 
206  points_tab = dal::dynamic_tas<ind_cv_ct, 8>();
207  convex_tab = dal::dynamic_tas<mesh_convex_structure, 8>();
208 
209  }
210 
211  void mesh_structure::stat(void) {
212  cout << "mesh structure with " << nb_convex() << " valid convex, "
213  << "for a total memory size of " << memsize() << " bytes.\n";
214  }
215 
216 
218  ind_set &s) const {
219  s.resize(0);
221 
222  for (size_type i = 0; i < points_tab[pt[0]].size(); ++i) {
223  size_type icv = points_tab[pt[0]][i];
224  if (icv != ic && is_convex_having_points(icv, short_type(pt.size()),
225  pt.begin())
226  && (convex_tab[ic].cstruct->dim()==convex_tab[icv].cstruct->dim()))
227  s.push_back(icv);
228  }
229  }
230 
231 
233  const std::vector<short_type> &ftab,
234  ind_set &s) const {
235  s.resize(0);
236  std::vector<size_type> ipts;
237 
238  switch (ftab.size()) {
239  case 0:
240  {
241  ind_cv_ct pt = ind_points_of_convex(ic);
242  ipts.resize(pt.size());
243  std::copy(pt.begin(), pt.end(), ipts.begin());
244  }
245  break;
246 
247  case 1:
248  {
250  ipts.resize(pt.size());
251  std::copy(pt.begin(), pt.end(), ipts.begin());
252  }
253  break;
254 
255  default:
256  {
257  const mesh_convex_structure &q = convex_tab[ic];
258  const convex_ind_ct &ind = q.cstruct->ind_common_points_of_faces(ftab);
259  if (ind.size() == 0) return neighbours_of_convex(ic, s);
260  ipts.resize(ind.size());
261  auto it = ind.cbegin();
262  for (size_type &ipt : ipts) ipt = q.pts[*it++];
263  }
264  break;
265  }
266 
267  if (ipts.size() == 0) {
268  return; // Should we return the all the neighbours ?
269  }
270 
271  auto ipt0 = ipts.cbegin();
272  auto ipt1 = ipt0 + 1;
273  short_type nbpts = short_type(ipts.size()-1);
274  for (size_type icv : points_tab[*ipt0]) {
275  if (icv != ic && (nbpts==0 || is_convex_having_points(icv, nbpts, ipt1)))
276  s.push_back(icv);
277  }
278  }
279 
280  void mesh_structure::neighbours_of_convex(size_type ic, ind_set &s) const {
281  s.resize(0);
282  unsigned nbf = nb_faces_of_convex(ic);
283  for (short_type iff = 0; iff < nbf; ++iff) {
285 
286  for (size_type i = 0; i < points_tab[pt[0]].size(); ++i) {
287  size_type icv = points_tab[pt[0]][i];
288  if (icv != ic && is_convex_having_points(icv, short_type(pt.size()),
289  pt.begin())
290  && (convex_tab[ic].cstruct->dim()==convex_tab[icv].cstruct->dim()))
291  if (std::find(s.begin(), s.end(), icv) == s.end())
292  s.push_back(icv);
293  }
294  }
295  }
296 
298  short_type iff) const {
300 
301  for (size_type i = 0; i < points_tab[pt[0]].size(); ++i) {
302  size_type icv = points_tab[pt[0]][i];
303  if (icv != ic && is_convex_having_points(icv, short_type(pt.size()),
304  pt.begin())
305  && (convex_tab[ic].cstruct->dim()==convex_tab[icv].cstruct->dim()))
306  return icv;
307  }
308  return size_type(-1);
309  }
310 
312  size_type neighbour_element = neighbour_of_convex(cv, f);
313  if (neighbour_element == size_type(-1)) return convex_face::invalid_face();
314  auto pcs = structure_of_convex(neighbour_element);
315  auto face_points = ind_points_of_face_of_convex(cv, f);
316  auto nNeighbourElementFaces = pcs->nb_faces();
317  for (short_type iff = 0; iff < nNeighbourElementFaces; ++iff) {
318  auto nPointsOnFace = pcs->nb_points_of_face(iff);
319  if (is_convex_face_having_points(neighbour_element, iff,
320  nPointsOnFace, face_points.begin()))
321  return {neighbour_element, iff};
322  }
323  GMM_ASSERT2(false, "failed to determine neighbouring face");
324  return convex_face::invalid_face();
325  }
326 
328  size_type ip) const {
329  const ind_cv_ct &ct = ind_points_of_convex(ic);
330  ind_cv_ct::const_iterator it = std::find(ct.begin(), ct.end(), ip);
331  return (it != ct.end()) ? (it - ct.begin()): size_type(-1);
332  }
333 
334  /* ********************************************************************* */
335  /* */
336  /* Gives the list of edges of a convex face. */
337  /* */
338  /* ********************************************************************* */
339 
340  void mesh_edge_list_convex(pconvex_structure cvs, std::vector<size_type> points_of_convex,
341  size_type cv_id, edge_list &el, bool merge_convex)
342  { // a tester ... optimisable.
343  dim_type n = cvs->dim();
344  short_type nbp = cvs->nb_points();
345  size_type ncv = merge_convex ? 0 : cv_id;
346 
347  if (nbp == n+1 && cvs == simplex_structure(n)) {
348  for (dim_type k = 0; k < n; ++k)
349  for (dim_type l = dim_type(k+1); l <= n; ++l)
350  el.add(edge_list_elt(points_of_convex[k],
351  points_of_convex[l], ncv));
352  }
353  else if (nbp == (size_type(1) << n)
354  && cvs == parallelepiped_structure(n)) {
355  for (size_type k = 0; k < (size_type(1) << n); ++k)
356  for (dim_type j = 0; j < n; ++j)
357  if ((k & (1 << j)) == 0)
358  el.add(edge_list_elt(points_of_convex[k],
359  points_of_convex[k | (1 << j)], ncv));
360  }
361  else if (nbp == 2 * n && cvs == prism_P1_structure(n)) {
362  for (dim_type k = 0; k < n - 1; ++k)
363  for (dim_type l = dim_type(k+1); l < n; ++l) {
364  el.add(edge_list_elt(points_of_convex[k],
365  points_of_convex[l], ncv));
366  el.add(edge_list_elt(points_of_convex[k+n],
367  points_of_convex[l+n], ncv));
368  }
369  for (dim_type k = 0; k < n; ++k)
370  el.add(edge_list_elt(points_of_convex[k],
371  points_of_convex[k+n], ncv));
372  }
373  else {
376  size_type ncs = 1;
377  cvstab[0] = cvs;
378  indpttab[0].resize(cvstab[0]->nb_points());
379  std::copy(points_of_convex.begin(),
380  points_of_convex.end(), indpttab[0].begin());
381 
382  /* pseudo recursive decomposition of the initial convex into
383  face of face of ... of face of convex , cvstab and indpttab
384  being the "stack" of convexes to treat
385  */
386  while (ncs != 0) {
387  ncs--;
388  cvs = cvstab[ncs];
389  std::vector< size_type > pts = indpttab[ncs];
390  if (cvs->dim() == 1) { // il faudrait ├ętendre aux autres cas classiques.
391 
392  for (size_type j = 1; j < cvs->nb_points(); ++j) {
393  //cerr << "ncs=" << ncs << "j=" << j << ", ajout de " << (indpttab[ncs])[j-1] << "," << (indpttab[ncs])[j] << endl;
394  el.add(edge_list_elt((indpttab[ncs])[j-1],(indpttab[ncs])[j],ncv));
395  }
396  }
397  else {
398  short_type nf = cvs->nb_faces();
399  //cerr << "ncs = " << ncs << ",cvs->dim=" << int(cvs->dim()) << ", nb_faces=" << nf << endl;
400  for (short_type f = 0; f < nf; ++f) {
401  cvstab[ncs+f] = (cvs->faces_structure())[f];
402  indpttab[ncs+f].resize(cvs->nb_points_of_face(f));
403  /*
404  cerr << " -> f=" << f << ", cvs->nb_points_of_face(f)=" <<
405  cvs->nb_points_of_face(f) << "=[";
406  for (size_type k = 0; k < cvs->nb_points_of_face(f); ++k)
407  cerr << (std::string(k==0?"":","))
408  << int(cvs->ind_points_of_face(f)[k]) << "{" << pts[cvs->ind_points_of_face(f)[k]] << "}";
409  cerr << "]" << endl;
410  */
411  for (size_type k = 0; k < cvs->nb_points_of_face(f); ++k)
412  (indpttab[ncs+f])[k] = pts[cvs->ind_points_of_face(f)[k]];
413  }
414  // cvstab[ncs] = cvstab[ncs + nf - 1];
415  //indpttab[ncs] = indpttab[ncs + nf - 1];
416  ncs += nf;
417  //cerr << "on empile les " << nf << " faces, -> ncs=" << ncs << endl;
418  }
419  }
420  }
421  }
422 
423 
424  void mesh_edge_list(const mesh_structure &m, edge_list &el,
425  bool merge_convex) {
426  std::vector<size_type> p;
427  for (dal::bv_visitor cv(m.convex_index()); !cv.finished(); ++cv) {
428  p.resize(m.nb_points_of_convex(cv));
429  std::copy(m.ind_points_of_convex(cv).begin(),
430  m.ind_points_of_convex(cv).end(), p.begin());
431  mesh_edge_list_convex(m.structure_of_convex(cv), p, cv,
432  el, merge_convex);
433  }
434  }
435 
436 
437  // static double enumerate_dof_time = 0;
439  std::vector<size_type> &cmk) {
440  const dal::bit_vector &cvlst = ms.convex_index();
441  cmk.reserve(cvlst.card()); cmk.resize(0);
442  if (ms.convex_index().card() == 0) return;
443  std::deque<int> pile;
444  std::vector<size_type> tab, connectivity(cvlst.last_true()+1),
445  temp(cvlst.last_true()+1);
446 
447  size_type cv = cvlst.first_true();
448 
449  std::fill(connectivity.begin(), connectivity.end(), size_type(-1));
450  // double t = dal::uclock_sec();
451 
452  /* count neighbours for each convex */
453  for (dal::bv_visitor j(cvlst); !j.finished(); ++j) {
454  const mesh_structure::ind_cv_ct &ct = ms.ind_points_of_convex(j);
455  mesh_structure::ind_cv_ct::const_iterator itp = ct.begin(),
456  itpe = ct.end();
457  size_type nei = 0;
458 
459  for (; itp != itpe; ++itp) {
460  size_type ip = *itp;
461  mesh_structure::ind_cv_ct::const_iterator
462  it = ms.convex_to_point(ip).begin(),
463  ite = ms.convex_to_point(ip).end();
464  for ( ; it != ite; ++it)
465  if (temp[*it] != j+1) { temp[*it] = j+1; nei++; }
466  }
467  connectivity[j] = nei-1;
468  if (nei < connectivity[cv]) cv = j;
469  }
470 
471  /* do the cuthill mckee */
472  while (cv != size_type(-1)) {
473  connectivity[cv] = size_type(-1);
474  cmk.push_back(cv);
475  size_type nbp = ms.nb_points_of_convex(cv);
476 
477  for (size_type i = 0; i < nbp; i++) {
478  size_type ip = ms.ind_points_of_convex(cv)[i];
479  mesh_structure::ind_cv_ct::const_iterator
480  it = ms.convex_to_point(ip).begin(),
481  ite = ms.convex_to_point(ip).end();
482  for ( ; it != ite; ++it)
483  if (connectivity[*it] != size_type(-1)) {
484  connectivity[*it] = size_type(-1);
485  pile.push_front(int(*it));
486  }
487  }
488 
489  if (pile.empty()) {
490  cv = std::min_element(connectivity.begin(), connectivity.end()) - connectivity.begin();
491  if (connectivity[cv] == size_type(-1)) break;
492  }
493  else { cv = pile.back(); pile.pop_back(); }
494  }
495 
496  // enumerate_dof_time += dal::uclock_sec() - t;
497  // cerr << "cuthill_mckee_on_convexes: " << enumerate_dof_time << " sec\n";
498  }
499 
500 } /* end of namespace bgeot. */
void swap_convex(size_type cv1, size_type cv2)
Exchange two convex IDs.
indexed array reference (given a container X, and a set of indexes I, this class provides a pseudo-co...
Definition: gmm_ref.h:289
short_type nb_points_of_convex(size_type ic) const
Return the number of points of convex ic.
ind_pt_face_ct ind_points_of_face_of_convex(size_type ic, short_type f) const
Return a container of the (global) point number for face f or convex ic.
void to_faces(dim_type n)
build a new mesh, such that its convexes are the faces of the convexes of the previous one ...
void to_edges()
build a new mesh, such that its convexes are the edges of the convexes of the previous one ...
void optimize_structure()
Reorder the convex IDs and point IDs, such that there is no hole in their numbering.
pconvex_structure prism_P1_structure(dim_type nc)
Give a pointer on the structures of a prism of dimension d.
convex_face adjacent_face(size_type ic, short_type f) const
Return a face of the neighbouring element that is adjacent to the given face.
pconvex_structure structure_of_convex(size_type ic) const
Return the pconvex_structure of the convex ic.
const ind_cv_ct & convex_to_point(size_type ip) const
Return a container of the convexes attached to point ip.
void swap_points(size_type i, size_type j)
Exchange two point IDs.
void clear()
erase the mesh
bool is_convex_face_having_points(size_type ic, short_type face_num, short_type nb, ITER pit) const
Return true if the face of the convex contains the given list of points.
void ind_points_to_point(size_type, ind_set &) const
Return a container of the points attached (via an edge) to point ip.
std::shared_ptr< const convex_structure > pconvex_structure
Pointer on a convex structure description.
size_type nb_convex() const
The total number of convexes in the mesh.
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
bool is_convex_having_points(size_type ic, short_type nb, ITER pit) const
Return true if the convex contains the listed points.
Dynamic Array.
Definition: dal_basic.h:47
size_type memsize(void) const
Gives the total memory occupied by the array.
Definition: dal_basic.h:273
size_type add_face_of_convex(size_type ic, short_type f)
Insert a new convex corresponding to face f of the convex ic.
Mesh structure definition.
size_type neighbour_of_convex(size_type ic, short_type f) const
Return a neighbour convex of a given convex face.
size_type size(void) const
Number of allocated elements.
Definition: dal_basic.h:219
size_type local_ind_of_convex_point(size_type ic, size_type ip) const
Return the "local" index for point ip of the mesh.
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:79
pconvex_structure parallelepiped_structure(dim_type nc, dim_type k)
Give a pointer on the structures of a parallelepiped of dimension d.
size_type add_convex(pconvex_structure cs, ITER ipts, bool *present=0)
Insert a new convex in the mesh_structure.
pconvex_structure simplex_structure(dim_type nc)
Give a pointer on the structures of a simplex of dimension d.
iterator end(void)
Iterator on the last + 1 element.
Definition: dal_basic.h:229
Mesh structure definition.
size_type ind_in_convex_of_point(size_type ic, size_type ip) const
Find the local index of the point of global index ip with respect to the convex cv.
Basic Geometric Tools.
bool is_convex_valid(size_type i)
Return true if i is in convex_index()
void cuthill_mckee_on_convexes(const bgeot::mesh_structure &ms, std::vector< size_type > &cmk)
Return the cuthill_mc_kee ordering on the convexes.
const ind_cv_ct & ind_points_of_convex(size_type ic) const
Return a container to the list of points attached to convex ic.
void sup_convex(size_type ic)
Remove the convex ic.
iterator begin(void)
Iterator on the first element.
Definition: dal_basic.h:225
short_type nb_faces_of_convex(size_type ic) const
Return the number of faces of convex ic.
void neighbours_of_convex(size_type ic, short_type f, ind_set &s) const
Return in s a list of neighbours of a given convex face.
void add_faces_of_convex(size_type ic)
Insert new convexes corresponding to the faces of the convex ic.
const dal::bit_vector & convex_index() const
Return the list of valid convex IDs.