GetFEM  5.4.2
getfem_mesh.cc
1 /*===========================================================================
2 
3  Copyright (C) 1999-2020 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 #include "getfem/bgeot_torus.h"
23 #include "getfem/dal_singleton.h"
24 #include "gmm/gmm_condition_number.h"
25 #include "getfem/getfem_mesh.h"
26 #include "getfem/getfem_integration.h"
27 
28 #if GETFEM_HAVE_METIS_OLD_API
29 extern "C" void METIS_PartGraphKway(int *, int *, int *, int *, int *, int *,
30  int *, int *, int *, int *, int *);
31 #elif GETFEM_HAVE_METIS
32 # include <metis.h>
33 #endif
34 
35 namespace getfem {
36 
37  gmm::uint64_type act_counter(void) {
38  static gmm::uint64_type c = gmm::uint64_type(1);
39  return ++c;
40  }
41 
42  void mesh::sup_convex_from_regions(size_type c) {
43  for (dal::bv_visitor i(valid_cvf_sets); !i.finished(); ++i)
44  cvf_sets[i].sup_all(c);
45  touch();
46  }
47 
48  void mesh::swap_convex_in_regions(size_type c1, size_type c2) {
49  for (dal::bv_visitor i(valid_cvf_sets); !i.finished(); ++i)
50  cvf_sets[i].swap_convex(c1, c2);
51  touch();
52  }
53 
54  void mesh::handle_region_refinement(size_type ic,
55  const std::vector<size_type> &icv,
56  bool refine) {
57 
58  bgeot::pgeometric_trans pgt = trans_of_convex(ic);
59  ind_set s;
60 
61  for (dal::bv_visitor ir(valid_cvf_sets); !ir.finished(); ++ir) {
62  mesh_region &r = cvf_sets[ir];
63 
64 
65  if (refine && r[ic].any()) {
66  if (r[ic][0])
67  for (size_type jc = 0; jc < icv.size(); ++jc)
68  r.add(icv[jc]);
69 
70  bgeot::geotrans_inv_convex giv;
71  for (short_type f = 0; f < pgt->structure()->nb_faces(); ++f) {
72  if (r[ic][f+1]) {
73 
74  for (size_type jc = 0; jc < icv.size(); ++jc) {
75  bgeot::pgeometric_trans pgtsub = trans_of_convex(icv[jc]);
76  for (short_type fsub = 0; fsub < pgtsub->structure()->nb_faces();
77  ++fsub) {
78  neighbors_of_convex(icv[jc], fsub, s);
79  ind_set::const_iterator it = s.begin(), ite = s.end();
80  bool found = false;
81  for (; it != ite; ++it)
82  if (std::find(icv.begin(), icv.end(), *it) != icv.end())
83  { found = true; break; }
84  if (found) continue;
85 
86  base_node pt, barycentre
87  =gmm::mean_value(pgtsub->convex_ref()->points_of_face(fsub));
88  pt = pgtsub->transform(barycentre, points_of_convex(icv[jc]));
89 
90  giv.init(points_of_convex(ic), pgt);
91  giv.invert(pt, barycentre);
92 
93 
94  if (gmm::abs(pgt->convex_ref()->is_in_face(f,barycentre)) < 0.001)
95  r.add(icv[jc], fsub);
96  }
97  }
98  }
99  }
100  }
101 
102  for (size_type jc = 0; jc < icv.size(); ++jc)
103  if (!refine && r[icv[jc]].any()) {
104  if (r[icv[jc]][0])
105  r.add(ic);
106  bgeot::geotrans_inv_convex giv;
107  bgeot::pgeometric_trans pgtsub = trans_of_convex(icv[jc]);
108  for (short_type fsub = 0; fsub < pgtsub->structure()->nb_faces();
109  ++fsub)
110  if (r[icv[jc]][fsub+1]) {
111  base_node pt, barycentre
112  = gmm::mean_value(pgtsub->convex_ref()->points_of_face(fsub));
113  pt = pgtsub->transform(barycentre, points_of_convex(icv[jc]));
114 
115  giv.init(points_of_convex(ic), pgt);
116  giv.invert(pt, barycentre);
117 
118  for (short_type f = 0; f < pgt->structure()->nb_faces(); ++f)
119  if (gmm::abs(pgt->convex_ref()->is_in_face(f,barycentre)) < 0.001)
120  { r.add(ic, f); break; }
121  }
122  }
123  }
124  }
125 
126  void mesh::init(void) {
127 #if GETFEM_PARA_LEVEL > 1
128  modified = true;
129 #endif
130  cuthill_mckee_uptodate = false;
131  }
132 
133  mesh::mesh(const std::string name) : name_(name) { init(); }
134 
135  mesh::mesh(const bgeot::basic_mesh &m, const std::string name)
136  : bgeot::basic_mesh(m), name_(name) { init(); }
137 
138  mesh::mesh(const mesh &m) : context_dependencies(),
139  std::enable_shared_from_this<getfem::mesh>()
140  { copy_from(m); }
141 
142  mesh &mesh::operator=(const mesh &m) {
143  clear_dependencies();
144  clear();
145  copy_from(m);
146  return *this;
147  }
148 
149 #if GETFEM_PARA_LEVEL > 1
150 
151  void mesh::compute_mpi_region() const {
152  int size, rank;
153  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
154  MPI_Comm_size(MPI_COMM_WORLD, &size);
155 
156  if (size < 2) {
157  mpi_region = mesh_region::all_convexes();
158  mpi_region.from_mesh(*this);
159  } else {
160  int ne = int(nb_convex());
161  std::vector<int> xadj(ne+1), adjncy, numelt(ne), npart(ne);
162  std::vector<int> indelt(nb_allocated_convex());
163 
164  double t_ref = MPI_Wtime();
165 
166  int j = 0, k = 0;
167  ind_set s;
168  for (dal::bv_visitor ic(convex_index()); !ic.finished(); ++ic, ++j) {
169  numelt[j] = ic;
170  indelt[ic] = j;
171  }
172  j = 0;
173  for (dal::bv_visitor ic(convex_index()); !ic.finished(); ++ic, ++j) {
174  xadj[j] = k;
175  neighbors_of_convex(ic, s);
176  for (ind_set::iterator it = s.begin();
177  it != s.end(); ++it) { adjncy.push_back(indelt[*it]); ++k; }
178  }
179  xadj[j] = k;
180 
181 #ifdef GETFEM_HAVE_METIS_OLD_API
182  int wgtflag = 0, numflag = 0, edgecut;
183  int options[5] = {0,0,0,0,0};
184  METIS_PartGraphKway(&ne, &(xadj[0]), &(adjncy[0]), 0, 0, &wgtflag,
185  &numflag, &size, options, &edgecut, &(npart[0]));
186 #else
187  int ncon = 1, edgecut;
188  int options[METIS_NOPTIONS] = { 0 };
189  METIS_SetDefaultOptions(options);
190  METIS_PartGraphKway(&ne, &ncon, &(xadj[0]), &(adjncy[0]), 0, 0, 0,
191  &size, 0, 0, options, &edgecut, &(npart[0]));
192 #endif
193 
194  for (size_type i = 0; i < size_type(ne); ++i)
195  if (npart[i] == rank) mpi_region.add(numelt[i]);
196 
197  if (MPI_IS_MASTER())
198  cout << "Partition time "<< MPI_Wtime()-t_ref << endl;
199  }
200  modified = false;
201  valid_sub_regions.clear();
202  }
203 
204  void mesh::compute_mpi_sub_region(size_type n) const {
205  if (valid_cvf_sets.is_in(n)) {
206  mpi_sub_region[n] = mesh_region::intersection(cvf_sets[n], mpi_region);
207  }
208  else
209  mpi_sub_region[n] = mesh_region();
210  valid_sub_regions.add(n);
211  }
212 
213  void mesh::intersect_with_mpi_region(mesh_region &rg) const {
214  if (rg.id() == mesh_region::all_convexes().id()) {
215  rg = get_mpi_region();
216  } else if (int(rg.id()) >= 0) {
217  rg = get_mpi_sub_region(rg.id());
218  } else
219  rg = mesh_region::intersection(rg, get_mpi_region());
220  }
221 #endif
222 
223  void mesh::optimize_structure(bool with_renumbering) {
224  pts.resort();
225  size_type i, j = nb_convex(), nbc = j;
226  for (i = 0; i < j; i++)
227  if (!convex_tab.index_valid(i))
228  swap_convex(i, convex_tab.ind_last());
229  if (pts.size())
230  for (i = 0, j = pts.size()-1;
231  i < j && j != ST_NIL; ++i, --j) {
232  while (i < j && j != ST_NIL && pts.index()[i]) ++i;
233  while (i < j && j != ST_NIL && !(pts.index()[j])) --j;
234  if (i < j && j != ST_NIL ) swap_points(i, j);
235  }
236  if (with_renumbering) { // Could be optimized no using only swap_convex
237  std::vector<size_type> cmk, iord(nbc), iordinv(nbc);
238  for (i = 0; i < nbc; ++i) iord[i] = iordinv[i] = i;
239 
240  bgeot::cuthill_mckee_on_convexes(*this, cmk);
241  for (i = 0; i < nbc; ++i) {
242  j = iordinv[cmk[i]];
243  if (i != j) {
244  swap_convex(i, j);
245  std::swap(iord[i], iord[j]);
246  std::swap(iordinv[iord[i]], iordinv[iord[j]]);
247  }
248  }
249  }
250  }
251 
252  void mesh::translation(const base_small_vector &V)
253  { pts.translation(V); touch(); }
254 
255  void mesh::transformation(const base_matrix &M) {
256  pts.transformation(M);
257  Bank_info = std::unique_ptr<Bank_info_struct>();
258  touch();
259  }
260 
261  void mesh::bounding_box(base_node& Pmin, base_node& Pmax) const {
262  bool is_first = true;
263  Pmin.clear(); Pmax.clear();
264  for (dal::bv_visitor i(pts.index()); !i.finished(); ++i) {
265  if (is_first) { Pmin = Pmax = pts[i]; is_first = false; }
266  else for (dim_type j=0; j < dim(); ++j) {
267  Pmin[j] = std::min(Pmin[j], pts[i][j]);
268  Pmax[j] = std::max(Pmax[j], pts[i][j]);
269  }
270  }
271  }
272 
273  void mesh::clear(void) {
274  mesh_structure::clear();
275  pts.clear();
276  gtab.clear(); trans_exists.clear();
277  cvf_sets.clear(); valid_cvf_sets.clear();
278  cvs_v_num.clear();
279  Bank_info = nullptr;
280  touch();
281  }
282 
283  size_type mesh::add_segment(size_type a, size_type b) {
284  size_type ipt[2]; ipt[0] = a; ipt[1] = b;
285  return add_convex(bgeot::simplex_geotrans(1, 1), &(ipt[0]));
286  }
287 
288  size_type mesh::add_triangle(size_type a,
289  size_type b, size_type c) {
290  size_type ipt[3]; ipt[0] = a; ipt[1] = b; ipt[2] = c;
291  return add_simplex(2, &(ipt[0]));
292  }
293 
294  size_type mesh::add_triangle_by_points
295  (const base_node &p1, const base_node &p2, const base_node &p3)
296  { return add_triangle(add_point(p1), add_point(p2), add_point(p3)); }
297 
298  size_type mesh::add_tetrahedron(size_type a, size_type b,
299  size_type c, size_type d) {
300  size_type ipt[4]; ipt[0] = a; ipt[1] = b; ipt[2] = c; ipt[3] = d;
301  return add_simplex(3, &(ipt[0]));
302  }
303 
304  size_type mesh::add_pyramid(size_type a, size_type b,
305  size_type c, size_type d, size_type e) {
306  size_type ipt[5] = {a, b, c, d, e};
307  return add_convex(bgeot::pyramid_QK_geotrans(1), &(ipt[0]));
308  }
309 
310  size_type mesh::add_tetrahedron_by_points
311  (const base_node &p1, const base_node &p2,
312  const base_node &p3, const base_node &p4) {
313  return add_tetrahedron(add_point(p1), add_point(p2),
314  add_point(p3), add_point(p4));
315  }
316 
317  void mesh::merge_convexes_from_mesh(const mesh &msource, size_type rg,
318  scalar_type tol) {
319 
320  size_type nbpt = points_index().last()+1;
321  GMM_ASSERT1(nbpt == nb_points(),
322  "Please call the optimize_structure() function before "
323  "merging elements from another mesh");
324  GMM_ASSERT1(rg == size_type(-1) || msource.region(rg).is_only_convexes(),
325  "The provided mesh region should only contain convexes");
326 
327  const dal::bit_vector &convexes = (rg == size_type(-1))
328  ? msource.convex_index()
329  : msource.region(rg).index();
330  std::vector<size_type> old2new(msource.points_index().last()+1, size_type(-1));
331  for (dal::bv_visitor cv(convexes); !cv.finished(); ++cv) {
332 
333  bgeot::pgeometric_trans pgt = msource.trans_of_convex(cv);
334  short_type nb = short_type(pgt->nb_points());
335  const ind_cv_ct &rct = msource.ind_points_of_convex(cv);
336  GMM_ASSERT1(nb == rct.size(), "Internal error");
337  std::vector<size_type> ind(nb);
338  for (short_type i = 0; i < nb; ++i) {
339  size_type old_pid = rct[i];
340  size_type new_pid = old2new[old_pid];
341  if (new_pid == size_type(-1)) {
342  size_type next_pid = points_index().last()+1;
343  base_node pt = msource.points()[old_pid];
344  new_pid = pts.add_node(pt, tol);
345  if (new_pid < next_pid && new_pid >= nbpt) {
346  // do not allow internal merging of nodes in the source mesh
347  new_pid = pts.add_node(pt, -1.);
348  GMM_ASSERT1(new_pid == next_pid, "Internal error");
349  }
350  old2new[old_pid] = new_pid;
351  }
352  ind[i] = new_pid;
353  }
354  add_convex(pgt, ind.begin());
355  }
356  }
357 
358  void mesh::sup_convex(size_type ic, bool sup_points) {
359  static std::vector<size_type> ipt;
360  if (sup_points) {
361  const ind_cv_ct &ct = ind_points_of_convex(ic);
362  ipt.assign(ct.begin(), ct.end());
363  }
364  bgeot::mesh_structure::sup_convex(ic);
365  if (sup_points)
366  for (const size_type &ip : ipt) sup_point(ip);
367  trans_exists.sup(ic);
368  sup_convex_from_regions(ic);
369  if (Bank_info.get()) Bank_sup_convex_from_green(ic);
370  touch();
371  }
372 
373  void mesh::swap_convex(size_type i, size_type j) {
374  if (i != j) {
375  bgeot::mesh_structure::swap_convex(i,j);
376  trans_exists.swap(i, j);
377  gtab.swap(i,j);
378  swap_convex_in_regions(i, j);
379  if (Bank_info.get()) Bank_swap_convex(i,j);
380  cvs_v_num[i] = cvs_v_num[j] = act_counter(); touch();
381  }
382  }
383 
384  base_small_vector mesh::normal_of_face_of_convex(size_type ic, short_type f,
385  const base_node &pt) const {
386  bgeot::pgeometric_trans pgt = trans_of_convex(ic);
387  base_matrix G(dim(),pgt->nb_points());
388  vectors_to_base_matrix(G, points_of_convex(ic));
389  bgeot::geotrans_interpolation_context c(trans_of_convex(ic), pt, G);
390  return bgeot::compute_normal(c, f);
391  }
392 
393 
394  base_small_vector mesh::normal_of_face_of_convex(size_type ic, short_type f,
395  size_type n) const {
396  bgeot::pgeometric_trans pgt = trans_of_convex(ic);
397  bgeot::pgeotrans_precomp pgp
398  = bgeot::geotrans_precomp(pgt, pgt->pgeometric_nodes(), 0);
399  base_matrix G;
400  vectors_to_base_matrix(G, points_of_convex(ic));
401  bgeot::geotrans_interpolation_context
402  c(pgp,pgt->structure()->ind_points_of_face(f)[n], G);
403  return bgeot::compute_normal(c, f);
404  }
405 
406  base_small_vector mesh::mean_normal_of_face_of_convex(size_type ic,
407  short_type f) const {
408  bgeot::pgeometric_trans pgt = trans_of_convex(ic);
409  base_matrix G; vectors_to_base_matrix(G,points_of_convex(ic));
410  base_small_vector ptmean(dim());
411  size_type nbpt = pgt->structure()->nb_points_of_face(f);
412  for (size_type i = 0; i < nbpt; ++i)
413  gmm::add(pgt->geometric_nodes()[pgt->structure()->ind_points_of_face(f)[i]], ptmean);
414  ptmean /= scalar_type(nbpt);
415  bgeot::geotrans_interpolation_context c(pgt,ptmean, G);
416  base_small_vector n = bgeot::compute_normal(c, f);
417  n /= gmm::vect_norm2(n);
418  return n;
419  }
420 
421  base_matrix mesh::local_basis_of_face_of_convex(size_type ic, short_type f,
422  const base_node &pt) const {
423  bgeot::pgeometric_trans pgt = trans_of_convex(ic);
424  base_matrix G(dim(),pgt->nb_points());
425  vectors_to_base_matrix(G,points_of_convex(ic));
426  bgeot::geotrans_interpolation_context c(trans_of_convex(ic), pt, G);
427  return bgeot::compute_local_basis(c, f);
428  }
429 
430  base_matrix mesh::local_basis_of_face_of_convex(size_type ic, short_type f,
431  size_type n) const {
432  bgeot::pgeometric_trans pgt = trans_of_convex(ic);
433  bgeot::pgeotrans_precomp pgp
434  = bgeot::geotrans_precomp(pgt, pgt->pgeometric_nodes(), 0);
435  base_matrix G(dim(),pgt->nb_points());
436  vectors_to_base_matrix(G,points_of_convex(ic));
437  bgeot::geotrans_interpolation_context
438  c(pgp,pgt->structure()->ind_points_of_face(f)[n], G);
439  return bgeot::compute_local_basis(c, f);
440  }
441 
442  scalar_type mesh::convex_area_estimate(size_type ic, size_type deg) const {
443  base_matrix G;
444  bgeot::vectors_to_base_matrix(G, points_of_convex(ic));
445  return getfem::convex_area_estimate
446  (trans_of_convex(ic), G, classical_approx_im(trans_of_convex(ic),
447  dim_type(deg)));
448  }
449 
450  scalar_type mesh::convex_quality_estimate(size_type ic) const {
451  base_matrix G;
452  bgeot::vectors_to_base_matrix(G, points_of_convex(ic));
453  auto pgt = trans_of_convex(ic);
454  if (auto pgt_torus = dynamic_cast<const bgeot::torus_geom_trans*>(pgt.get())) {
455  pgt = pgt_torus->get_original_transformation();
456  G.resize(pgt->dim(), G.ncols());
457  }
458  return getfem::convex_quality_estimate(pgt, G);
459  }
460 
461  scalar_type mesh::convex_radius_estimate(size_type ic) const {
462  base_matrix G;
463  bgeot::vectors_to_base_matrix(G, points_of_convex(ic));
464  return getfem::convex_radius_estimate(trans_of_convex(ic), G);
465  }
466 
467  scalar_type mesh::minimal_convex_radius_estimate() const {
468  if (convex_index().empty()) return 1;
469  scalar_type r = convex_radius_estimate(convex_index().first_true());
470  for (dal::bv_visitor cv(convex_index()); !cv.finished(); ++cv) {
471  r = std::min(r, convex_radius_estimate(cv));
472  }
473  return r;
474  }
475 
476  scalar_type mesh::maximal_convex_radius_estimate() const {
477  if (convex_index().empty()) return 1;
478  scalar_type r = convex_radius_estimate(convex_index().first_true());
479  for (dal::bv_visitor cv(convex_index()); !cv.finished(); ++cv) {
480  r = std::max(r, convex_radius_estimate(cv));
481  }
482  return r;
483  }
484 
485  void mesh::set_name(const std::string& name){name_=name;}
486 
487  void mesh::copy_from(const mesh& m) {
488  clear();
489  set_name(m.name_);
490  bgeot::basic_mesh::operator=(m);
491  for (const auto &kv : m.cvf_sets) {
492  if (kv.second.get_parent_mesh() != 0)
493  cvf_sets[kv.first].set_parent_mesh(this);
494  cvf_sets[kv.first] = kv.second;
495  }
496  valid_cvf_sets = m.valid_cvf_sets;
497  cvs_v_num.clear();
498  gmm::uint64_type d = act_counter();
499  for (dal::bv_visitor i(convex_index()); !i.finished(); ++i)
500  cvs_v_num[i] = d;
501  Bank_info = std::unique_ptr<Bank_info_struct>();
502  if (m.Bank_info.get())
503  Bank_info = std::make_unique<Bank_info_struct>(*(m.Bank_info));
504  }
505 
506  struct mesh_convex_structure_loc {
507  bgeot::pgeometric_trans cstruct;
508  std::vector<size_type> pts;
509  };
510 
511  void mesh::read_from_file(std::istream &ist) {
512  gmm::stream_standard_locale sl(ist);
513  dal::bit_vector npt;
514  dal::dynamic_array<double> tmpv;
515  std::string tmp;
516  bool te = false, please_get = true;
517 
518  ist.precision(16);
519  clear();
520  ist.seekg(0);ist.clear();
521  bgeot::read_until(ist, "BEGIN POINTS LIST");
522 
523  while (!te) {
524  if (please_get) bgeot::get_token(ist, tmp); else please_get = true;
525 
526  if (!bgeot::casecmp(tmp, "END"))
527  { te = true; }
528  else if (!bgeot::casecmp(tmp, "POINT")) {
529  bgeot::get_token(ist, tmp);
530  if (!bgeot::casecmp(tmp, "COUNT")) {
531  bgeot::get_token(ist, tmp); // Ignored. Used in some applications
532  } else {
533  size_type ip = atoi(tmp.c_str());
534  dim_type d = 0;
535  GMM_ASSERT1(!npt.is_in(ip),
536  "Two points with the same index. loading aborted.");
537  npt.add(ip);
538  bgeot::get_token(ist, tmp);
539  while (isdigit(tmp[0]) || tmp[0] == '-' || tmp[0] == '+'
540  || tmp[0] == '.')
541  { tmpv[d++] = atof(tmp.c_str()); bgeot::get_token(ist, tmp); }
542  please_get = false;
543  base_node v(d);
544  for (size_type i = 0; i < d; i++) v[i] = tmpv[i];
545  size_type ipl = add_point(v);
546  if (ip != ipl) {
547  GMM_ASSERT1(!npt.is_in(ipl), "Two points [#" << ip << " and #"
548  << ipl << "] with the same coords "<< v
549  << ". loading aborted.");
550  swap_points(ip, ipl);
551  }
552  }
553  } else if (tmp.size()) {
554  GMM_ASSERT1(false, "Syntax error in file, at token '" << tmp
555  << "', pos=" << std::streamoff(ist.tellg()));
556  } else if (ist.eof()) {
557  GMM_ASSERT1(false, "Unexpected end of stream while reading mesh");
558  }
559  }
560 
561  bool tend = false;
562  dal::dynamic_array<mesh_convex_structure_loc> cv;
563  dal::bit_vector ncv;
564 
565  ist.seekg(0);
566  if (!bgeot::read_until(ist, "BEGIN MESH STRUCTURE DESCRIPTION"))
567  GMM_ASSERT1(false, "This seems not to be a mesh file");
568 
569  while (!tend) {
570  tend = !bgeot::get_token(ist, tmp);
571  if (!bgeot::casecmp(tmp, "END"))
572  { tend = true; }
573  else if (!bgeot::casecmp(tmp, "CONVEX")) {
574  size_type ic;
575  bgeot::get_token(ist, tmp);
576  if (!bgeot::casecmp(tmp, "COUNT")) {
577  bgeot::get_token(ist, tmp); // Ignored. Used in some applications
578  } else {
579  ic = gmm::abs(atoi(tmp.c_str()));
580  GMM_ASSERT1(!ncv.is_in(ic),
581  "Negative or repeated index, loading aborted.");
582  ncv.add(ic);
583 
584  int rgt = bgeot::get_token(ist, tmp);
585  if (rgt != 3) { // for backward compatibility with version 1.7
586  char c; ist.get(c);
587  while (!isspace(c)) { tmp.push_back(c); ist.get(c); }
588  }
589 
590  bgeot::pgeometric_trans pgt = bgeot::geometric_trans_descriptor(tmp);
591  size_type nb = pgt->nb_points();
592 
593  cv[ic].cstruct = pgt;
594  cv[ic].pts.resize(nb);
595  for (size_type i = 0; i < nb; i++) {
596  bgeot::get_token(ist, tmp);
597  cv[ic].pts[i] = gmm::abs(atoi(tmp.c_str()));
598  }
599  }
600  }
601  else if (tmp.size()) {
602  GMM_ASSERT1(false, "Syntax error reading a mesh file "
603  " at pos " << std::streamoff(ist.tellg())
604  << "(expecting 'CONVEX' or 'END', found '"<< tmp << "')");
605  } else if (ist.eof()) {
606  GMM_ASSERT1(false, "Unexpected end of stream "
607  << "(missing BEGIN MESH/END MESH ?)");
608  }
609  }
610  ist >> bgeot::skip("MESH STRUCTURE DESCRIPTION");
611 
612  for (dal::bv_visitor ic(ncv); !ic.finished(); ++ic) {
613  size_type i = add_convex(cv[ic].cstruct, cv[ic].pts.begin());
614  if (i != ic) swap_convex(i, ic);
615  }
616 
617  tend = false;
618  while (!tend) {
619  tend = !bgeot::get_token(ist, tmp);
620  // bool error = false;
621  if (bgeot::casecmp(tmp, "BEGIN")==0) {
622  bgeot::get_token(ist, tmp);
623  if (bgeot::casecmp(tmp, "BOUNDARY")==0 ||
624  bgeot::casecmp(tmp, "REGION")==0) {
625  bgeot::get_token(ist, tmp);
626  size_type bnum = atoi(tmp.c_str());
627  bgeot::get_token(ist, tmp);
628  while (true) {
629  if (bgeot::casecmp(tmp, "END")!=0) {
630  size_type ic = atoi(tmp.c_str());
631  bgeot::get_token(ist, tmp);
632  if (tmp[0] == '/') {
633  bgeot::get_token(ist, tmp);
634  if (!bgeot::casecmp(tmp, "END")) break;
635  int f = atoi(tmp.c_str());
636  region(bnum).add(ic, short_type(f));
637  bgeot::get_token(ist, tmp);
638  }
639  else {
640  region(bnum).add(ic);
641  if (!bgeot::casecmp(tmp, "END")) break;
642  }
643  } else break;
644  }
645  bgeot::get_token(ist, tmp);
646  bgeot::get_token(ist, tmp);
647  } else tend = true;
648  /*else GMM_ASSERT1(false, "Syntax error in file at token '"
649  << tmp << "' [pos=" << std::streamoff(ist.tellg())
650  << "]");*/
651  } else tend=true;
652  }
653  }
654 
655  void mesh::read_from_file(const std::string &name) {
656  std::ifstream o(name.c_str());
657  GMM_ASSERT1(o, "Mesh file '" << name << "' does not exist");
658  read_from_file(o);
659  o.close();
660  }
661 
662  template<class ITER>
663  void write_tab_to_file_(std::ostream &ost, const ITER& b_, const ITER& e)
664  { for (ITER b(b_) ; b != e; ++b) ost << " " << *b; }
665 
666  template<class ITER>
667  static void write_convex_to_file_(const mesh &ms,
668  std::ostream &ost,
669  ITER b, ITER e) {
670  for ( ; b != e; ++b) {
671  size_type i = b.index();
672  ost << " CONVEX " << i << " \'"
673  << bgeot::name_of_geometric_trans(ms.trans_of_convex(i)).c_str()
674  << "\' ";
675  write_tab_to_file_(ost, ms.ind_points_of_convex(i).begin(),
676  ms.ind_points_of_convex(i).end() );
677  ost << '\n';
678  }
679  }
680 
681  template<class ITER> static void write_point_to_file_(std::ostream &ost,
682  ITER b, ITER e)
683  { for ( ; b != e; ++b) ost << " " << *b; ost << '\n'; }
684 
685  void mesh::write_to_file(std::ostream &ost) const {
686  ost.precision(16);
687  gmm::stream_standard_locale sl(ost);
688  ost << '\n' << "BEGIN POINTS LIST" << '\n' << '\n';
689  ost << " POINT COUNT " << points().index().last_true()+1 << '\n';
690  for (size_type i = 0; i < points_tab.size(); ++i)
691  if (is_point_valid(i) ) {
692  ost << " POINT " << i;
693  write_point_to_file_(ost, pts[i].begin(), pts[i].end());
694  }
695  ost << '\n' << "END POINTS LIST" << '\n' << '\n' << '\n';
696 
697  ost << '\n' << "BEGIN MESH STRUCTURE DESCRIPTION" << '\n' << '\n';
698  ost << " CONVEX COUNT " << convex_index().last_true()+1 << '\n';
699  write_convex_to_file_(*this, ost, convex_tab.tas_begin(),
700  convex_tab.tas_end());
701  ost << '\n' << "END MESH STRUCTURE DESCRIPTION" << '\n';
702 
703  for (dal::bv_visitor bnum(valid_cvf_sets); !bnum.finished(); ++bnum) {
704  ost << "BEGIN REGION " << bnum << "\n" << region(bnum) << "\n"
705  << "END REGION " << bnum << "\n";
706  }
707  }
708 
709  void mesh::write_to_file(const std::string &name) const {
710  std::ofstream o(name.c_str());
711  GMM_ASSERT1(o, "impossible to write to file '" << name << "'");
712  o << "% GETFEM MESH FILE " << '\n';
713  o << "% GETFEM VERSION " << GETFEM_VERSION << '\n' << '\n' << '\n';
714  write_to_file(o);
715  o.close();
716  }
717 
718  size_type mesh::memsize(void) const {
719  return bgeot::mesh_structure::memsize() - sizeof(bgeot::mesh_structure)
720  + pts.memsize() + (pts.index().last_true()+1)*dim()*sizeof(scalar_type)
721  + sizeof(mesh) + trans_exists.memsize() + gtab.memsize()
722  + valid_cvf_sets.card()*sizeof(mesh_region) + valid_cvf_sets.memsize();
723  }
724 
725  struct equilateral_to_GT_PK_grad_aux : public std::vector<base_matrix> {};
726  static const base_matrix &equilateral_to_GT_PK_grad(dim_type N) {
727  std::vector<base_matrix>
728  &pbm = dal::singleton<equilateral_to_GT_PK_grad_aux >::instance();
729  if (N > pbm.size()) pbm.resize(N);
730  if (pbm[N-1].empty()) {
731  bgeot::pgeometric_trans pgt = bgeot::simplex_geotrans(N,1);
732  base_matrix Gr(N,N);
733  base_matrix G(N,N+1);
734  vectors_to_base_matrix
735  (G, bgeot::equilateral_simplex_of_reference(N)->points());
736  gmm::mult(G, bgeot::geotrans_precomp
737  (pgt, pgt->pgeometric_nodes(), 0)->grad(0), Gr);
738  gmm::lu_inverse(Gr);
739  pbm[N-1].swap(Gr);
740  }
741  return pbm[N-1];
742  }
743 
744 
745  scalar_type convex_area_estimate(bgeot::pgeometric_trans pgt,
746  const base_matrix& G,
747  pintegration_method pi) {
748  double area(0);
749  static bgeot::pgeometric_trans pgt_old = 0;
750  static bgeot::pgeotrans_precomp pgp = 0;
751  static pintegration_method pim_old = 0;
752  papprox_integration pai = get_approx_im_or_fail(pi);
753  if (pgt_old != pgt || pim_old != pi) {
754  pgt_old = pgt;
755  pim_old = pi;
756  pgp = bgeot::geotrans_precomp
757  (pgt, pai->pintegration_points(), pi);
758  }
759  bgeot::geotrans_interpolation_context gic(pgp, 0, G);
760  for (size_type i = 0; i < pai->nb_points_on_convex(); ++i) {
761  gic.set_ii(i);
762  area += pai->coeff(i) * gic.J();
763  }
764  return area;
765  }
766 
767  /* TODO : use the geotrans from an "equilateral" reference element to
768  the real element
769  check if the sign of the determinants does change
770  (=> very very bad quality of the convex)
771  */
772  scalar_type convex_quality_estimate(bgeot::pgeometric_trans pgt,
773  const base_matrix& G) {
774  static bgeot::pgeometric_trans pgt_old = 0;
775  static bgeot::pgeotrans_precomp pgp = 0;
776  if (pgt_old != pgt) {
777  pgt_old=pgt;
778  pgp=bgeot::geotrans_precomp(pgt, pgt->pgeometric_nodes(), 0);
779  }
780 
781  size_type n = (pgt->is_linear()) ? 1 : pgt->nb_points();
782  scalar_type q = 1;
783  dim_type N = dim_type(G.nrows()), P = pgt->structure()->dim();
784  base_matrix K(N,P);
785  for (size_type ip=0; ip < n; ++ip) {
786  gmm::mult(G, pgp->grad(ip), K);
787  /* TODO : this is an ugly fix for simplexes only.. there should be
788  a transformation of any pgt to the equivalent equilateral pgt
789  (for prisms etc) */
790  if (bgeot::basic_structure(pgt->structure())
791  == bgeot::simplex_structure(P))
792  gmm::mult(base_matrix(K),equilateral_to_GT_PK_grad(P),K);
793  q = std::max(q, gmm::condition_number(K));
794  }
795  return 1./q;
796  }
797 
798  scalar_type convex_radius_estimate(bgeot::pgeometric_trans pgt,
799  const base_matrix& G) {
800  static bgeot::pgeometric_trans pgt_old = 0;
801  static bgeot::pgeotrans_precomp pgp = 0;
802  if (pgt_old != pgt) {
803  pgt_old=pgt;
804  pgp=bgeot::geotrans_precomp(pgt, pgt->pgeometric_nodes(), 0);
805  }
806  size_type N = G.nrows();
807  size_type n = (pgt->is_linear()) ? 1 : pgt->nb_points();
808  scalar_type q = 0;
809  base_matrix K(pgp->grad(0).ncols(),G.nrows());
810  for (size_type ip=0; ip < n; ++ip) {
811  gmm::mult(gmm::transposed(pgp->grad(ip)), gmm::transposed(G), K);
812  scalar_type emax, emin; gmm::condition_number(K,emax,emin);
813  q = std::max(q, emax);
814  }
815  return q * sqrt(scalar_type(N)) / scalar_type(2);
816  }
817 
818  /* extract faces of convexes which are not shared
819  + convexes whose dimension is smaller that m.dim()
820  */
821  void
822  outer_faces_of_mesh(const getfem::mesh &m, const dal::bit_vector& cvlst,
823  convex_face_ct& flist) {
824  for (dal::bv_visitor ic(cvlst); !ic.finished(); ++ic) {
825  if (m.structure_of_convex(ic)->dim() == m.dim()) {
826  for (short_type f = 0; f < m.structure_of_convex(ic)->nb_faces(); f++) {
827  if (m.neighbor_of_convex(ic,f) == size_type(-1)) {
828  flist.push_back(convex_face(ic,f));
829  }
830  }
831  } else {
832  flist.push_back(convex_face(ic, short_type(-1)));
833  }
834  }
835  }
836 
837  void outer_faces_of_mesh(const mesh &m,
838  const mesh_region &cvlst,
839  mesh_region &flist) {
840  cvlst.error_if_not_convexes();
841  for (mr_visitor i(cvlst); !i.finished(); ++i) {
842  if (m.structure_of_convex(i.cv())->dim() == m.dim()) {
843  for (short_type f = 0; f < m.structure_of_convex(i.cv())->nb_faces();
844  f++) {
845  size_type cv2 = m.neighbor_of_convex(i.cv(), f);
846  if (cv2 == size_type(-1) || !cvlst.is_in(cv2)) {
847  flist.add(i.cv(),f);
848  }
849  }
850  }
851  else flist.add(i.cv());
852  }
853  }
854 
855  /* Select all the faces of the given mesh region (counted twice if they
856  are shared by two neighbor elements)
857  */
858  mesh_region all_faces_of_mesh(const mesh &m, const mesh_region &mr) {
859  mesh_region mrr;
860  mr.from_mesh(m);
861  mr.error_if_not_convexes();
862 
863  for (mr_visitor i(mr); !i.finished(); ++i) {
864  size_type cv1 = i.cv();
865  short_type nbf = m.structure_of_convex(i.cv())->nb_faces();
866  for (short_type f = 0; f < nbf; ++f)
867  mrr.add(cv1, f);
868  }
869  return mrr;
870  }
871 
872  /* Select all the faces sharing at least two element of the given mesh
873  region. Each face is represented only once and is arbitrarily chosen
874  between the two neighbor elements. Try to minimize the number of
875  elements.
876  */
877  mesh_region inner_faces_of_mesh(const mesh &m, const mesh_region &mr) {
878  mesh_region mrr;
879  mr.from_mesh(m);
880  mr.error_if_not_convexes();
881  dal::bit_vector visited;
882  bgeot::mesh_structure::ind_set neighbors;
883 
884  for (mr_visitor i(mr); !i.finished(); ++i) {
885  size_type cv1 = i.cv();
886  short_type nbf = m.structure_of_convex(i.cv())->nb_faces();
887  bool neighbor_visited = false;
888  for (short_type f = 0; f < nbf; ++f) {
889  neighbors.resize(0); m.neighbors_of_convex(cv1, f, neighbors);
890  for (size_type j = 0; j < neighbors.size(); ++j)
891  if (visited.is_in(neighbors[j]))
892  { neighbor_visited = true; break; }
893  }
894  if (!neighbor_visited) {
895  for (short_type f = 0; f < nbf; ++f) {
896  size_type cv2 = m.neighbor_of_convex(cv1, f);
897  if (cv2 != size_type(-1) && mr.is_in(cv2)) mrr.add(cv1,f);
898  }
899  visited.add(cv1);
900  }
901  }
902 
903  for (mr_visitor i(mr); !i.finished(); ++i) {
904  size_type cv1 = i.cv();
905  if (!(visited.is_in(cv1))) {
906  short_type nbf = m.structure_of_convex(i.cv())->nb_faces();
907  for (short_type f = 0; f < nbf; ++f) {
908  neighbors.resize(0); m.neighbors_of_convex(cv1, f, neighbors);
909  bool ok = false;
910  for (size_type j = 0; j < neighbors.size(); ++j) {
911  if (visited.is_in(neighbors[j])) { ok = false; break; }
912  if (mr.is_in(neighbors[j])) { ok = true; }
913  }
914  if (ok) { mrr.add(cv1,f); }
915  }
916  visited.add(cv1);
917  }
918  }
919  return mrr;
920  }
921 
922 
923  mesh_region select_faces_of_normal(const mesh &m, const mesh_region &mr,
924  const base_small_vector &V,
925  scalar_type angle) {
926  mesh_region mrr;
927  scalar_type threshold = gmm::vect_norm2(V)*cos(angle);
928  for (getfem::mr_visitor i(mr); !i.finished(); ++i)
929  if (i.f() != short_type(-1)) {
930  base_node un = m.mean_normal_of_face_of_convex(i.cv(), i.f());
931  if (gmm::vect_sp(V, un) - threshold >= -1E-8)
932  mrr.add(i.cv(), i.f());
933  }
934  return mrr;
935  }
936 
937  mesh_region select_faces_in_box(const mesh &m, const mesh_region &mr,
938  const base_node &pt1, const base_node &pt2) {
939  mesh_region mrr;
940  size_type N = m.dim();
941  GMM_ASSERT1(pt1.size() == N && pt2.size() == N, "Wrong dimensions");
942  for (getfem::mr_visitor i(mr, m); !i.finished(); ++i)
943  if (i.f() != short_type(-1)) {
944  bgeot::mesh_structure::ind_pt_face_ct pt
945  = m.ind_points_of_face_of_convex(i.cv(), i.f());
946 
947  bool is_in = true;
948  for (bgeot::mesh_structure::ind_pt_face_ct::iterator it = pt.begin();
949  it != pt.end(); ++it) {
950  for (size_type j = 0; j < N; ++j)
951  if (m.points()[*it][j] < pt1[j] || m.points()[*it][j] > pt2[j])
952  { is_in = false; break; }
953  if (!is_in) break;
954  }
955  if (is_in) mrr.add(i.cv(), i.f());
956  }
957  return mrr;
958  }
959 
960  mesh_region select_faces_in_ball(const mesh &m, const mesh_region &mr,
961  const base_node &center, scalar_type radius) {
962  mesh_region mrr;
963  size_type N = m.dim();
964  GMM_ASSERT1(center.size() == N, "Wrong dimensions");
965  for (getfem::mr_visitor i(mr, m); !i.finished(); ++i)
966  if (i.f() != short_type(-1)) {
967  bgeot::mesh_structure::ind_pt_face_ct pt
968  = m.ind_points_of_face_of_convex(i.cv(), i.f());
969 
970  bool is_in = true;
971  for (bgeot::mesh_structure::ind_pt_face_ct::iterator it = pt.begin();
972  it != pt.end(); ++it) {
973  scalar_type checked_radius = scalar_type(0.0);
974  for (size_type j = 0; j < N; ++j)
975  checked_radius += pow(m.points()[*it][j] - center[j], 2);
976  checked_radius = std::sqrt(checked_radius);
977  if (checked_radius > radius) { is_in = false; break; }
978  }
979  if (is_in) mrr.add(i.cv(), i.f());
980  }
981  return mrr;
982  }
983 
984  mesh_region select_convexes_in_box(const mesh &m, const mesh_region &mr,
985  const base_node &pt1, const base_node &pt2) {
986  mesh_region mrr;
987  size_type N = m.dim();
988  GMM_ASSERT1(pt1.size() == N && pt2.size() == N, "Wrong dimensions");
989  for (getfem::mr_visitor i(mr, m); !i.finished(); ++i)
990  if (i.f() == short_type(-1)) {
991  bgeot::mesh_structure::ind_cv_ct pt = m.ind_points_of_convex(i.cv());
992 
993  bool is_in = true;
994  for (bgeot::mesh_structure::ind_cv_ct::iterator it = pt.begin();
995  it != pt.end(); ++it) {
996  for (size_type j = 0; j < N; ++j)
997  if (m.points()[*it][j] < pt1[j] || m.points()[*it][j] > pt2[j])
998  { is_in = false; break; }
999  if (!is_in) break;
1000  }
1001  if (is_in) mrr.add(i.cv());
1002  }
1003  return mrr;
1004  }
1005 
1006  void extrude(const mesh& in, mesh& out, size_type nb_layers, short_type degree) {
1007  dim_type dim = in.dim();
1008  base_node pt(dim+1);
1009  out.clear();
1010  size_type nbpt = in.points().index().last()+1;
1011  GMM_ASSERT1(nbpt == in.points().index().card(),
1012  "please call the optimize_structure() method before using "
1013  "the mesh as a base for prismatic mesh");
1014  size_type nb_layers_total = nb_layers * degree;
1015  for (const base_node &inpt : in.points()) {
1016  std::copy(inpt.begin(), inpt.end(), pt.begin());
1017  pt[dim] = 0.0;
1018  for (size_type j = 0; j <= nb_layers_total;
1019  ++j, pt[dim] += scalar_type(1) / scalar_type(nb_layers_total))
1020  out.add_point(pt);
1021  }
1022 
1023  std::vector<size_type> tab;
1024  for (dal::bv_visitor cv(in.convex_index()); !cv.finished(); ++cv) {
1025  size_type nbp = in.nb_points_of_convex(cv);
1026  tab.resize((degree+1)*nbp);
1027  for (size_type j = 0; j < nb_layers; ++j) {
1028  for (short_type d = 0; d < degree+1; ++d)
1029  for (size_type k = 0; k < nbp; ++k)
1030  tab[k+nbp*d] = (nb_layers_total+1)*in.ind_points_of_convex(cv)[k] + j*degree + d;
1031  bgeot::pgeometric_trans pgt =
1032  bgeot::product_geotrans(in.trans_of_convex(cv),
1033  bgeot::simplex_geotrans(1,degree));
1034  out.add_convex(pgt, tab.begin());
1035  }
1036  }
1037  }
1038 }
1039 
1040 
1041 //
1042 // Bank refinement
1043 //
1044 
1045 
1046 namespace bgeot {
1047  size_type refinement_simplexe_tab(size_type n, size_type **tab);
1048 }
1049 
1050 namespace getfem {
1051 
1052  bool mesh::edge::operator <(const edge &e) const {
1053  if (i0 < e.i0) return true;
1054  if (i0 > e.i0) return false;
1055  if (i1 < e.i1) return true;
1056  if (i1 > e.i1) return false;
1057  if (i2 < e.i2) return true;
1058  return false;
1059  }
1060 
1061  void mesh::Bank_sup_convex_from_green(size_type i) {
1062  if (Bank_info.get() && Bank_info->is_green_simplex.is_in(i)) {
1063  size_type igs = Bank_info->num_green_simplex[i];
1064  green_simplex &gs = Bank_info->green_simplices[igs];
1065  for (size_type j = 0; j < gs.sub_simplices.size(); ++j) {
1066  Bank_info->num_green_simplex.erase(gs.sub_simplices[j]);
1067  Bank_info->is_green_simplex.sup(gs.sub_simplices[j]);
1068  }
1069  Bank_info->green_simplices.sup(igs);
1070  }
1071  }
1072 
1073  void mesh::Bank_swap_convex(size_type i, size_type j) {
1074  if (Bank_info.get()) {
1075  Bank_info->is_green_simplex.swap(i, j);
1076  std::map<size_type, size_type>::iterator
1077  iti = Bank_info->num_green_simplex.find(i);
1078  std::map<size_type, size_type>::iterator
1079  ite = Bank_info->num_green_simplex.end();
1080  std::map<size_type, size_type>::iterator
1081  itj = Bank_info->num_green_simplex.find(j);
1082  size_type numi(0), numj(0);
1083  if (iti != ite)
1084  { numi = iti->second; Bank_info->num_green_simplex.erase(i); }
1085  if (itj != ite)
1086  { numj = itj->second; Bank_info->num_green_simplex.erase(j); }
1087  if (iti != ite) {
1088  Bank_info->num_green_simplex[j] = numi;
1089  green_simplex &gs = Bank_info->green_simplices[numi];
1090  for (size_type k = 0; k < gs.sub_simplices.size(); ++k)
1091  if (gs.sub_simplices[k] == i) gs.sub_simplices[k] = j;
1092  else if (gs.sub_simplices[k] == j) gs.sub_simplices[k] = i;
1093  }
1094  if (itj != ite) {
1095  Bank_info->num_green_simplex[i] = numj;
1096  if (iti == ite || numi != numj) {
1097  green_simplex &gs = Bank_info->green_simplices[numj];
1098  for (size_type k = 0; k < gs.sub_simplices.size(); ++k)
1099  if (gs.sub_simplices[k] == i) gs.sub_simplices[k] = j;
1100  else if (gs.sub_simplices[k] == j) gs.sub_simplices[k] = i;
1101  }
1102  }
1103  }
1104  }
1105 
1106  void mesh::Bank_build_first_mesh(mesh &m, size_type n) {
1107  bgeot::pconvex_ref pcr = bgeot::simplex_of_reference(dim_type(n), 2);
1108  m.clear();
1109  for (size_type ip = 0; ip < pcr->nb_points(); ++ip)
1110  m.add_point(pcr->points()[ip]);
1111  size_type *tab;
1112  size_type nbc = bgeot::refinement_simplexe_tab(n, &tab);
1113  for (size_type ic = 0; ic < nbc; ++ic, tab+=(n+1))
1114  m.add_simplex(dim_type(n), tab);
1115  }
1116 
1117  struct mesh_cache_for_Bank_basic_refine_convex : public mesh {};
1118 
1119  void mesh::Bank_basic_refine_convex(size_type i) {
1120  bgeot::pgeometric_trans pgt = trans_of_convex(i);
1121  size_type n = pgt->basic_structure()->dim();
1122 
1123  static bgeot::pgeometric_trans pgt1 = 0;
1124 
1125  mesh &mesh2 = dal::singleton<mesh_cache_for_Bank_basic_refine_convex>::instance();
1126 
1127  static bgeot::pstored_point_tab pspt = 0;
1128  static bgeot::pgeotrans_precomp pgp = 0;
1129  static std::vector<size_type> ipt, ipt2, icl;
1130 
1131  if (pgt != pgt1) {
1132  pgt1 = pgt;
1133  mesh mesh1;
1134  Bank_build_first_mesh(mesh1, n);
1135 
1136  mesh2.clear();
1137  ipt.resize(pgt->nb_points());
1138  for (size_type ic = 0; ic < mesh1.nb_convex(); ++ic) {
1139  assert(mesh1.convex_index().is_in(ic));
1140  bgeot::pgeometric_trans pgt2 = mesh1.trans_of_convex(ic);
1141  for (size_type ip = 0; ip < pgt->nb_points(); ++ip)
1142  ipt[ip] =
1143  mesh2.add_point(pgt2->transform(pgt->geometric_nodes()[ip],
1144  mesh1.points_of_convex(ic)));
1145  mesh2.add_convex(pgt, &ipt[0]);
1146  }
1147 
1148  // if (pspt) dal::del_stored_object(pspt);
1149  pspt = bgeot::store_point_tab(mesh2.points());
1150  pgp = bgeot::geotrans_precomp(pgt, pspt, 0);
1151  }
1152 
1153  base_node pt(n);
1154  ipt.resize(pspt->size());
1155  for (size_type ip = 0; ip < pspt->size(); ++ip) {
1156  pgp->transform(points_of_convex(i), ip, pt);
1157  ipt[ip] = add_point(pt);
1158  }
1159 
1160  ipt2.resize(pgt->nb_points()); icl.resize(0);
1161  for (size_type ic = 0; ic < mesh2.nb_convex(); ++ic) {
1162  for (size_type j = 0; j < pgt->nb_points(); ++j)
1163  ipt2[j] = ipt[mesh2.ind_points_of_convex(ic)[j]];
1164  icl.push_back(add_convex(pgt, ipt2.begin()));
1165  }
1166  handle_region_refinement(i, icl, true);
1167  sup_convex(i, true);
1168  }
1169 
1170  void mesh::Bank_convex_with_edge(size_type i1, size_type i2,
1171  std::vector<size_type> &ipt) {
1172  ipt.resize(0);
1173  for (size_type k = 0; k < points_tab[i1].size(); ++k) {
1174  size_type cv = points_tab[i1][k], found = 0;
1175  const std::vector<size_type> &loc_ind = trans_of_convex(cv)->vertices();
1176  for (size_type i = 0; i < loc_ind.size(); ++i) {
1177  size_type ipp = convex_tab[cv].pts[loc_ind[i]];
1178  if (ipp == i1) ++found;
1179  if (ipp == i2) ++found;
1180  }
1181  GMM_ASSERT1(found <= 2, "Invalid convex with repeated nodes ");
1182  if (found == 2) ipt.push_back(cv);
1183  }
1184  }
1185 
1186  bool mesh::Bank_is_convex_having_points(size_type cv,
1187  const std::vector<size_type> &ipt) {
1188  size_type found = 0;
1189  const std::vector<size_type> &loc_ind = trans_of_convex(cv)->vertices();
1190  for (size_type i = 0; i < loc_ind.size(); ++i)
1191  if (std::find(ipt.begin(), ipt.end(),
1192  convex_tab[cv].pts[loc_ind[i]]) != ipt.end()) ++found;
1193  return (found >= ipt.size());
1194  }
1195 
1196  void mesh::Bank_refine_normal_convex(size_type i) {
1197  bgeot::pgeometric_trans pgt = trans_of_convex(i);
1198  GMM_ASSERT1(pgt->basic_structure() == bgeot::simplex_structure(pgt->dim()),
1199  "Sorry, refinement is only working with simplices.");
1200 
1201  const std::vector<size_type> &loc_ind = pgt->vertices();
1202  for (size_type ip1 = 0; ip1 < loc_ind.size(); ++ip1)
1203  for (size_type ip2 = ip1+1; ip2 < loc_ind.size(); ++ip2)
1204  Bank_info->edges.insert(edge(ind_points_of_convex(i)[loc_ind[ip1]],
1205  ind_points_of_convex(i)[loc_ind[ip2]]));
1206  Bank_basic_refine_convex(i);
1207  }
1208 
1209  size_type mesh::Bank_test_and_refine_convex(size_type i,
1210  dal::bit_vector &b, bool ref) {
1211  if (Bank_info->is_green_simplex[i]) {
1212  size_type igs = Bank_info->num_green_simplex[i];
1213  green_simplex &gs = Bank_info->green_simplices[igs];
1214 
1215  size_type icc = add_convex_by_points(gs.pgt, gs.cv.points().begin());
1216  handle_region_refinement(icc, gs.sub_simplices, false);
1217  for (size_type ic = 0; ic < gs.sub_simplices.size(); ++ic) {
1218  sup_convex(gs.sub_simplices[ic], true);
1219  b.sup(gs.sub_simplices[ic]);
1220  }
1221  if (!ref)
1222  for (size_type ip = 0; ip < gs.ipt_loc.size(); ++ip)
1223  for (size_type jp = ip+1; jp < gs.ipt_loc.size(); ++jp)
1224  Bank_info->edges.insert
1225  (edge(ind_points_of_convex(icc)[gs.ipt_loc[ip]],
1226  ind_points_of_convex(icc)[gs.ipt_loc[jp]]));
1227  Bank_sup_convex_from_green(i);
1228  if (ref) { Bank_refine_normal_convex(icc); return size_type(-1); }
1229  else return icc;
1230  }
1231  else if (ref) Bank_refine_normal_convex(i);
1232  return size_type(-1);
1233  }
1234 
1235  struct mesh_cache_for_Bank_build_green_simplexes : public mesh {};
1236 
1237  void mesh::Bank_build_green_simplexes(size_type ic,
1238  std::vector<size_type> &ipt) {
1239  GMM_ASSERT1(convex_index().is_in(ic), "Internal error");
1240  size_type igs = Bank_info->green_simplices.add(green_simplex());
1241  green_simplex &gs(Bank_info->green_simplices[igs]);
1242  std::vector<base_node> pt_tab(nb_points_of_convex(ic));
1243  ref_mesh_pt_ct ptab = points_of_convex(ic);
1244  pt_tab.assign(ptab.begin(), ptab.end());
1245  gs.cv = bgeot::convex<base_node>(structure_of_convex(ic), pt_tab);
1246 
1247  bgeot::pgeometric_trans pgt = gs.pgt = trans_of_convex(ic);
1248 
1249  size_type d = ipt.size() - 1, n = structure_of_convex(ic)->dim();
1250  static size_type d0 = 0;
1251  static bgeot::pstored_point_tab pspt1 = 0;
1252  mesh &mesh1
1253  = dal::singleton<mesh_cache_for_Bank_build_green_simplexes>::instance();
1254  if (d0 != d) {
1255  d0 = d;
1256  Bank_build_first_mesh(mesh1, d);
1257  pspt1 = bgeot::store_point_tab(mesh1.points());
1258  }
1259 
1260  const std::vector<size_type> &loc_ind = pgt->vertices();
1261  const bgeot::mesh_structure::ind_cv_ct &ct = ind_points_of_convex(ic);
1262 
1263  gs.ipt_loc.resize(ipt.size());
1264  std::vector<size_type> ipt_other;
1265  size_type nb_found = 0;
1266  for (size_type ip = 0; ip < loc_ind.size(); ++ip) {
1267  bool found = false;
1268  for (size_type i = 0; i < ipt.size(); ++i)
1269  if (ct[loc_ind[ip]] == ipt[i])
1270  { gs.ipt_loc[i] = ip; found = true; ++nb_found; break; }
1271  if (!found) ipt_other.push_back(ip);
1272  }
1273  assert(nb_found == ipt.size());
1274 
1275  mesh mesh2;
1276  for (size_type ip = 0; ip < loc_ind.size(); ++ip)
1277  mesh2.add_point(pgt->geometric_nodes()[loc_ind[ip]]);
1278  size_type ic1 = mesh2.add_simplex(dim_type(d), gs.ipt_loc.begin());
1279  bgeot::pgeometric_trans pgt1 = mesh2.trans_of_convex(ic1);
1280  for (size_type i = 0; i < ipt.size(); ++i)
1281  gs.ipt_loc[i] = loc_ind[gs.ipt_loc[i]];
1282 
1283  bgeot::pgeotrans_precomp pgp = bgeot::geotrans_precomp(pgt1, pspt1, 0);
1284 
1285  std::vector<size_type> ipt1(pspt1->size());
1286  base_node pt(n);
1287  for (size_type i = 0; i < pspt1->size(); ++i) {
1288  pgp->transform(mesh2.points_of_convex(ic1), i, pt);
1289  ipt1[i] = mesh2.add_point(pt);
1290  }
1291  mesh2.sup_convex(ic1);
1292 
1293  std::vector<size_type> ipt2(n+1);
1294  for (size_type i = 0; i < mesh1.nb_convex(); ++i) {
1295  for (size_type j = 0; j <= d; ++j)
1296  ipt2[j] = ipt1[mesh1.ind_points_of_convex(i)[j]];
1297  for (size_type j = d+1; j <= n; ++j)
1298  ipt2[j] = ipt_other[j-d-1];
1299  mesh2.add_simplex(dim_type(n), ipt2.begin());
1300  }
1301 
1302  mesh mesh3;
1303  ipt1.resize(pgt->nb_points());
1304  for (dal::bv_visitor i(mesh2.convex_index()); !i.finished(); ++i) {
1305  bgeot::pgeometric_trans pgt2 = mesh2.trans_of_convex(i);
1306  for (size_type ip = 0; ip < pgt->nb_points(); ++ip)
1307  ipt1[ip] =
1308  mesh3.add_point(pgt2->transform(pgt->geometric_nodes()[ip],
1309  mesh2.points_of_convex(i)));
1310  mesh3.add_convex(pgt, ipt1.begin());
1311  }
1312 
1313 
1314  bgeot::pstored_point_tab pspt3 = bgeot::store_point_tab(mesh3.points());
1315  pgp = bgeot::geotrans_precomp(pgt, pspt3, 0);
1316 
1317  ipt1.resize(pspt3->size());
1318  for (size_type ip = 0; ip < pspt3->size(); ++ip) {
1319  pgp->transform(points_of_convex(ic), ip, pt);
1320  ipt1[ip] = add_point(pt);
1321  }
1322  // dal::del_stored_object(pspt3);
1323 
1324  ipt2.resize(pgt->nb_points());
1325  for (size_type icc = 0; icc < mesh3.nb_convex(); ++icc) {
1326  for (size_type j = 0; j < pgt->nb_points(); ++j)
1327  ipt2[j] = ipt1[mesh3.ind_points_of_convex(icc)[j]];
1328  size_type i = add_convex(pgt, ipt2.begin());
1329  gs.sub_simplices.push_back(i);
1330  Bank_info->is_green_simplex.add(i);
1331  Bank_info->num_green_simplex[i] = igs;
1332  }
1333 
1334  for (size_type ip1 = 0; ip1 < ipt.size(); ++ip1)
1335  for (size_type ip2 = ip1+1; ip2 < ipt.size(); ++ip2)
1336  Bank_info->edges.insert(edge(ipt[ip1], ipt[ip2]));
1337 
1338  handle_region_refinement(ic, gs.sub_simplices, true);
1339  sup_convex(ic, true);
1340  }
1341 
1342  void mesh::Bank_refine(dal::bit_vector b) {
1343  if (!(Bank_info.get())) Bank_info = std::make_unique<Bank_info_struct>();
1344 
1345  b &= convex_index();
1346  if (b.card() == 0) return;
1347 
1348  Bank_info->edges.clear();
1349  while (b.card() > 0)
1350  Bank_test_and_refine_convex(b.take_first(), b);
1351 
1352  std::vector<size_type> ipt;
1353  edge_set marked_convexes;
1354  while (Bank_info->edges.size()) {
1355  marked_convexes.clear();
1356  b = convex_index();
1357  edge_set::const_iterator it = Bank_info->edges.begin();
1358  edge_set::const_iterator ite = Bank_info->edges.end(), it2=it;
1359  assert(!Bank_info->edges.empty());
1360  for (; it != ite; it = it2) {
1361  if (it2 != ite) ++it2;
1362  Bank_convex_with_edge(it->i1, it->i2, ipt);
1363  if (ipt.size() == 0) ; // Bank_info->edges.erase(it);
1364  else for (size_type ic = 0; ic < ipt.size(); ++ic)
1365  marked_convexes.insert(edge(ipt[ic], it->i1, it->i2));
1366  }
1367 
1368  it = marked_convexes.begin(); ite = marked_convexes.end();
1369  if (it == ite) break;
1370 
1371  while (it != ite) {
1372  size_type ic = it->i0;
1373  ipt.resize(0);
1374  while (it != ite && it->i0 == ic) {
1375  bool found1 = false, found2 = false;
1376  for (size_type j = 0; j < ipt.size(); ++j) {
1377  if (ipt[j] == it->i1) found1 = true;
1378  if (ipt[j] == it->i2) found2 = true;
1379  }
1380  if (!found1) ipt.push_back(it->i1);
1381  if (!found2) ipt.push_back(it->i2);
1382  ++it;
1383  }
1384  if (b.is_in(ic)) {
1385  if (ipt.size() > structure_of_convex(ic)->dim())
1386  Bank_test_and_refine_convex(ic, b);
1387  else if (Bank_info->is_green_simplex[ic]) {
1388  size_type icc = Bank_test_and_refine_convex(ic, b, false);
1389  if (!Bank_is_convex_having_points(icc, ipt)) {
1390  Bank_test_and_refine_convex(icc, b);
1391  }
1392  }
1393  else Bank_build_green_simplexes(ic, ipt);
1394  }
1395  }
1396  }
1397  Bank_info->edges.clear();
1398  }
1399 
1400  struct dummy_mesh_ {
1401  mesh m;
1402  dummy_mesh_() : m() {}
1403  };
1404 
1405  const mesh &dummy_mesh()
1406  { return dal::singleton<dummy_mesh_>::instance().m; }
1407 
1408 } /* end of namespace getfem. */
dal::singleton::instance
static T & instance()
Instance from the current thread.
Definition: dal_singleton.h:165
getfem::mesh::local_basis_of_face_of_convex
base_matrix local_basis_of_face_of_convex(size_type ic, short_type f, const base_node &pt) const
Return a local basis for the convex face.
Definition: getfem_mesh.cc:421
getfem::mesh::clear
void clear()
Erase the mesh.
Definition: getfem_mesh.cc:273
getfem::mesh::write_to_file
void write_to_file(const std::string &name) const
Write the mesh to a file.
Definition: getfem_mesh.cc:709
bgeot::compute_normal
base_small_vector compute_normal(const geotrans_interpolation_context &c, size_type face)
norm of returned vector is the ratio between the face surface on the real element and the face surfac...
Definition: bgeot_geometric_trans.cc:1082
bgeot::name_of_geometric_trans
std::string name_of_geometric_trans(pgeometric_trans p)
Get the string name of a geometric transformation.
Definition: bgeot_geometric_trans.cc:1173
bgeot_torus.h
Provides mesh of torus.
getfem::mesh::read_from_file
void read_from_file(const std::string &name)
Load the mesh from a file.
Definition: getfem_mesh.cc:655
bgeot::geotrans_inv_convex
does the inversion of the geometric transformation for a given convex
Definition: bgeot_geotrans_inv.h:64
getfem::mesh::swap_convex
void swap_convex(size_type i, size_type j)
Swap the indexes of the convex of indexes i and j in the whole structure.
Definition: getfem_mesh.cc:373
getfem::mesh::convex_area_estimate
scalar_type convex_area_estimate(size_type ic, size_type degree=2) const
Return an estimate of the convex area.
Definition: getfem_mesh.cc:442
getfem::mesh::region
const mesh_region region(size_type id) const
Return the region of index 'id'.
Definition: getfem_mesh.h:421
getfem::mesh::add_simplex
size_type add_simplex(dim_type di, ITER ipts)
Add a simplex to the mesh, given its dimension and point numbers.
Definition: getfem_mesh.h:252
bgeot::get_token
int get_token(std::istream &ist, std::string &st, bool ignore_cr, bool to_up, bool read_un_pm, int *linenb)
Very simple lexical analysis of general interest for reading small languages with a "MATLAB like" syn...
Definition: bgeot_ftool.cc:50
getfem::mesh::add_tetrahedron
size_type add_tetrahedron(size_type a, size_type b, size_type c, size_type d)
Add a tetrahedron to the mesh, given the point id of its vertices.
Definition: getfem_mesh.cc:298
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::convex_area_estimate
scalar_type APIDECL convex_area_estimate(bgeot::pgeometric_trans pgt, const base_matrix &pts, pintegration_method pim)
rough estimate of the convex area.
Definition: getfem_mesh.cc:745
bgeot::geotrans_interpolation_context::J
scalar_type J() const
get the Jacobian of the geometric trans (taken at point xref() )
Definition: bgeot_geometric_trans.h:453
getfem::inner_faces_of_mesh
mesh_region APIDECL inner_faces_of_mesh(const mesh &m, const mesh_region &mr=mesh_region::all_convexes())
Select all the faces sharing at least two element of the given mesh region.
Definition: getfem_mesh.cc:877
getfem_integration.h
Integration methods (exact and approximated) on convexes.
getfem::select_faces_in_ball
mesh_region APIDECL select_faces_in_ball(const mesh &m, const mesh_region &mr, const base_node &center, scalar_type radius)
Select in the region mr the faces of the mesh m lying entirely in the ball delimated by pt1 and radiu...
Definition: getfem_mesh.cc:960
bgeot::geotrans_interpolation_context
the geotrans_interpolation_context structure is passed as the argument of geometric transformation in...
Definition: bgeot_geometric_trans.h:411
bgeot::convex< base_node >
bgeot::geometric_trans_descriptor
pgeometric_trans geometric_trans_descriptor(std::string name)
Get the geometric transformation from its string name.
Definition: bgeot_geometric_trans.cc:1163
getfem::mesh::add_convex
size_type add_convex(bgeot::pgeometric_trans pgt, ITER ipts)
Add a convex to the mesh.
Definition: getfem_mesh.h:226
getfem::select_faces_in_box
mesh_region APIDECL select_faces_in_box(const mesh &m, const mesh_region &mr, const base_node &pt1, const base_node &pt2)
Select in the region mr the faces of the mesh m lying entirely in the box delimated by pt1 and pt2.
Definition: getfem_mesh.cc:937
getfem::convex_radius_estimate
scalar_type APIDECL convex_radius_estimate(bgeot::pgeometric_trans pgt, const base_matrix &pts)
rough estimate of the radius of the convex using the largest eigenvalue of the jacobian of the geomet...
Definition: getfem_mesh.cc:798
getfem::mesh::convex_quality_estimate
scalar_type convex_quality_estimate(size_type ic) const
Return an estimate of the convex quality (0 <= Q <= 1).
Definition: getfem_mesh.cc:450
dal::dynamic_array::begin
iterator begin(void)
Iterator on the first element.
Definition: dal_basic.h:225
getfem::mesh::normal_of_face_of_convex
base_small_vector normal_of_face_of_convex(size_type ic, short_type f, const base_node &pt) const
Return the normal of the given convex face, evaluated at the point pt.
Definition: getfem_mesh.cc:384
getfem::mesh::add_triangle_by_points
size_type add_triangle_by_points(const base_node &p1, const base_node &p2, const base_node &p3)
Add a triangle to the mesh, given the coordinates of its vertices.
Definition: getfem_mesh.cc:295
bgeot::mesh_structure::is_point_valid
bool is_point_valid(size_type i) const
Return true if the point i is used by at least one convex.
Definition: bgeot_mesh_structure.h:104
bgeot::mesh_structure::ind_points_of_convex
const ind_set & ind_points_of_convex(size_type ic) const
Return a container to the list of points attached to convex ic.
Definition: bgeot_mesh_structure.h:107
bgeot::basic_structure
pconvex_structure basic_structure(pconvex_structure cv)
Original structure (if concerned)
Definition: bgeot_convex_structure.h:172
getfem_mesh.h
Define a getfem::getfem_mesh object.
getfem::mesh::merge_convexes_from_mesh
void merge_convexes_from_mesh(const mesh &m, size_type rg=size_type(-1), scalar_type tol=scalar_type(0))
Merge all convexes from a another mesh, possibly restricted to a mesh region.
Definition: getfem_mesh.cc:317
getfem::mesh_region::all_convexes
static mesh_region all_convexes()
provide a default value for the mesh_region parameters of assembly procedures etc.
Definition: getfem_mesh_region.h:142
getfem::mesh::sup_point
void sup_point(size_type i)
Delete the point of index i from the mesh if it is not linked to a convex.
Definition: getfem_mesh.h:200
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
bgeot::simplex_structure
pconvex_structure simplex_structure(dim_type nc)
Give a pointer on the structures of a simplex of dimension d.
Definition: bgeot_convex_structure.cc:148
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
bgeot::mesh_structure::neighbor_of_convex
size_type neighbor_of_convex(size_type ic, short_type f) const
Return a neighbor convex of a given convex face.
Definition: bgeot_mesh_structure.cc:297
gmm::tab_ref_index_ref_iterator_
iterator over a gmm::tab_ref_index_ref<ITER,ITER_INDEX>
Definition: gmm_ref.h:226
getfem::mesh_region::intersection
static mesh_region intersection(const mesh_region &a, const mesh_region &b)
return the intersection of two mesh regions
Definition: getfem_mesh_region.cc:391
dal::dynamic_array::clear
void clear(void)
Clear and desallocate all the elements.
Definition: dal_basic.h:298
getfem::mesh::copy_from
void copy_from(const mesh &m)
Clone a mesh.
Definition: getfem_mesh.cc:487
getfem::mesh::convex_radius_estimate
virtual scalar_type convex_radius_estimate(size_type ic) const
Return an estimate of the convex largest dimension.
Definition: getfem_mesh.cc:461
getfem::mesh::sup_convex_from_regions
void sup_convex_from_regions(size_type cv)
Remove all references to a convex from all regions stored in the mesh.
Definition: getfem_mesh.cc:42
getfem::mesh_region::visitor
"iterator" class for regions.
Definition: getfem_mesh_region.h:237
getfem::mesh_region::index
const dal::bit_vector & index() const
Index of the region convexes, or the convexes from the partition on the current thread.
Definition: getfem_mesh_region.cc:237
getfem::mesh::points_index
const dal::bit_vector & points_index() const
Return the points index.
Definition: getfem_mesh.h:196
dal::dynamic_array
Dynamic Array.
Definition: dal_basic.h:47
bgeot::mesh_structure::ind_points_of_face_of_convex
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.
Definition: bgeot_mesh_structure.cc:173
bgeot::geotrans_interpolation_context::set_ii
void set_ii(size_type ii__)
change the current point (assuming a geotrans_precomp_ is used)
Definition: bgeot_geometric_trans.h:463
bgeot::mesh_structure::optimize_structure
void optimize_structure()
Reorder the convex IDs and point IDs, such that there is no hole in their numbering.
Definition: bgeot_mesh_structure.cc:191
getfem::mesh::bounding_box
void bounding_box(base_node &Pmin, base_node &Pmax) const
Return the bounding box [Pmin - Pmax] of the mesh.
Definition: getfem_mesh.cc:261
getfem::mesh_region::from_mesh
const mesh_region & from_mesh(const mesh &m) const
For regions which have been built with just a number 'id', from_mesh(m) sets the current region to 'm...
Definition: getfem_mesh_region.cc:72
bgeot::mesh_structure::nb_convex
size_type nb_convex() const
The total number of convexes in the mesh.
Definition: bgeot_mesh_structure.h:96
getfem::mesh::mesh
mesh(const std::string name="")
Constructor.
Definition: getfem_mesh.cc:133
bgeot::mesh_structure::swap_convex
void swap_convex(size_type cv1, size_type cv2)
Exchange two convex IDs.
Definition: bgeot_mesh_structure.cc:69
getfem::mesh::transformation
void transformation(const base_matrix &)
apply the given matrix transformation to each mesh node.
Definition: getfem_mesh.cc:255
getfem::mesh::translation
void translation(const base_small_vector &)
Apply the given translation to each mesh node.
Definition: getfem_mesh.cc:252
getfem::mesh::add_convex_by_points
size_type add_convex_by_points(bgeot::pgeometric_trans pgt, ITER ipts, const scalar_type tol=scalar_type(0))
Add a convex to the mesh, given a geometric transformation and a list of point coordinates.
Definition: getfem_mesh.h:558
getfem::mesh::sup_convex
void sup_convex(size_type ic, bool sup_points=false)
Delete the convex of index ic from the mesh.
Definition: getfem_mesh.cc:358
getfem::mesh::Bank_refine
void Bank_refine(dal::bit_vector)
Use the Bank strategy to refine some convexes.
Definition: getfem_mesh.cc:1342
getfem::mesh::add_segment
size_type add_segment(size_type a, size_type b)
Add a segment to the mesh, given the point id of its vertices.
Definition: getfem_mesh.cc:283
bgeot::node_tab::add_node
size_type add_node(const base_node &pt, const scalar_type radius=0, bool remove_duplicated_nodes=true)
Add a point to the array or use an existing point, located within a distance smaller than radius.
Definition: bgeot_node_tab.cc:96
getfem::mesh::add_tetrahedron_by_points
size_type add_tetrahedron_by_points(const base_node &p1, const base_node &p2, const base_node &p3, const base_node &p4)
Add a tetrahedron to the mesh, given the coordinates of its vertices.
Definition: getfem_mesh.cc:311
bgeot::mesh_structure
Mesh structure definition.
Definition: bgeot_mesh_structure.h:73
getfem::classical_approx_im
pintegration_method classical_approx_im(bgeot::pgeometric_trans pgt, dim_type degree)
try to find an approximate integration method for the geometric transformation pgt which is able to i...
Definition: getfem_integration.cc:1280
getfem::mesh::nb_points
size_type nb_points() const
Give the number of geometrical nodes in the mesh.
Definition: getfem_mesh.h:194
bgeot
Basic Geometric Tools.
Definition: bgeot_convex_ref.cc:27
bgeot::mesh_structure::nb_allocated_convex
size_type nb_allocated_convex() const
The number of convex indexes from 0 to the index of the last convex.
Definition: bgeot_mesh_structure.h:98
getfem::convex_quality_estimate
scalar_type APIDECL convex_quality_estimate(bgeot::pgeometric_trans pgt, const base_matrix &pts)
rough estimate of the maximum value of the condition number of the jacobian of the geometric transfor...
Definition: getfem_mesh.cc:772
getfem::mesh::add_triangle
size_type add_triangle(size_type a, size_type b, size_type c)
Add a triangle to the mesh, given the point id of its vertices.
Definition: getfem_mesh.cc:288
getfem::mesh
Describe a mesh (collection of convexes (elements) and points).
Definition: getfem_mesh.h:95
bgeot::torus_geom_trans
An adaptor that adapts a two dimensional geometric_trans to include radial dimension.
Definition: bgeot_torus.h:48
getfem::mesh::add_pyramid
size_type add_pyramid(size_type a, size_type b, size_type c, size_type d, size_type e)
Add a pyramid to the mesh, given the point id of its vertices.
Definition: getfem_mesh.cc:304
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:186
getfem::mesh_region
structure used to hold a set of convexes and/or convex faces.
Definition: getfem_mesh_region.h:55
dal::dynamic_array::size
size_type size(void) const
Number of allocated elements.
Definition: dal_basic.h:219
bgeot::geotrans_inv_convex::invert
bool invert(const base_node &n, base_node &n_ref, scalar_type IN_EPS=1e-12, bool project_into_element=false)
given the node on the real element, returns the node on the reference element (even if it is outside ...
Definition: bgeot_geotrans_inv.cc:61
bgeot::mesh_structure::sup_convex
void sup_convex(size_type ic)
Remove the convex ic.
Definition: bgeot_mesh_structure.cc:100
getfem::mesh_region::is_only_convexes
bool is_only_convexes() const
Return true if the region do not contain any convex face.
Definition: getfem_mesh_region.cc:345
getfem::all_faces_of_mesh
mesh_region APIDECL all_faces_of_mesh(const mesh &m, const mesh_region &mr=mesh_region::all_convexes())
Select all the faces of the given mesh region.
Definition: getfem_mesh.cc:858
getfem::outer_faces_of_mesh
void APIDECL outer_faces_of_mesh(const mesh &m, const dal::bit_vector &cvlst, convex_face_ct &flist)
returns a list of "exterior" faces of a mesh (i.e.
Definition: getfem_mesh.cc:822
bgeot::simplex_of_reference
pconvex_ref simplex_of_reference(dim_type nc, short_type K)
returns a simplex of reference of dimension nc and degree k
Definition: bgeot_convex_ref.cc:340
getfem::select_faces_of_normal
mesh_region APIDECL select_faces_of_normal(const mesh &m, const mesh_region &mr, const base_small_vector &V, scalar_type angle)
Select in the region mr the faces of the mesh m with their unit outward vector having a maximal angle...
Definition: getfem_mesh.cc:923
gmm_condition_number.h
computation of the condition number of dense matrices.
bgeot::equilateral_simplex_of_reference
pconvex_ref equilateral_simplex_of_reference(dim_type nc)
equilateral simplex (degree 1).
Definition: bgeot_convex_ref.cc:831
getfem::mesh::swap_points
void swap_points(size_type i, size_type j)
Swap the indexes of points of index i and j in the whole structure.
Definition: getfem_mesh.h:202
gmm::tab_ref_index_ref
indexed array reference (given a container X, and a set of indexes I, this class provides a pseudo-co...
Definition: gmm_ref.h:289
bgeot::mesh_structure::nb_points_of_convex
short_type nb_points_of_convex(size_type ic) const
Return the number of points of convex ic.
Definition: bgeot_mesh_structure.h:115
getfem::extrude
void APIDECL extrude(const mesh &in, mesh &out, size_type nb_layers, short_type degree=short_type(1))
build a N+1 dimensions mesh from a N-dimensions mesh by extrusion.
Definition: getfem_mesh.cc:1006
getfem::mesh::minimal_convex_radius_estimate
scalar_type minimal_convex_radius_estimate() const
Return an estimate of the convex smallest dimension.
Definition: getfem_mesh.cc:467
bgeot::mesh_structure::structure_of_convex
pconvex_structure structure_of_convex(size_type ic) const
Return the pconvex_structure of the convex ic.
Definition: bgeot_mesh_structure.h:112
getfem::mesh::maximal_convex_radius_estimate
scalar_type maximal_convex_radius_estimate() const
Return an estimate of the convex largest dimension.
Definition: getfem_mesh.cc:476
bgeot::cuthill_mckee_on_convexes
void cuthill_mckee_on_convexes(const bgeot::mesh_structure &ms, std::vector< size_type > &cmk)
Return the cuthill_mc_kee ordering on the convexes.
Definition: bgeot_mesh_structure.cc:438
bgeot::compute_local_basis
base_matrix compute_local_basis(const geotrans_interpolation_context &c, size_type face)
return the local basis (i.e.
Definition: bgeot_geometric_trans.cc:1094
bgeot::mesh_structure::neighbors_of_convex
void neighbors_of_convex(size_type ic, short_type f, ind_set &s) const
Return in s a list of neighbors of a given convex face.
Definition: bgeot_mesh_structure.cc:217