GetFEM  5.4.2
getfem_mesh_region.cc
1 /*===========================================================================
2 
3  Copyright (C) 2005-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 
23 #include "getfem/getfem_mesh.h"
24 #include "getfem/getfem_omp.h"
25 
26 namespace getfem {
27 
28  using face_bitset = mesh_region::face_bitset;
29 
30  mesh_region::mesh_region(const mesh_region &other)
31  : p(std::make_shared<impl>()), id_(size_type(-2)), parent_mesh(0) {
32  this->operator=(other);
33  }
34 
35  mesh_region::mesh_region()
36  : p(std::make_shared<impl>()), id_(size_type(-2)), type_(size_type(-1)),
37  partitioning_allowed{true}, parent_mesh(nullptr){
38  if (me_is_multithreaded_now()) prohibit_partitioning();
39  mark_region_changed();
40  }
41 
42  mesh_region::mesh_region(size_type id__) : id_(id__), type_(size_type(-1)),
43  partitioning_allowed{true}, parent_mesh(nullptr){
44  mark_region_changed();
45  }
46 
47  mesh_region::mesh_region(mesh& m, size_type id__, size_type type) :
48  p(std::make_shared<impl>()), id_(id__), type_(type),
49  partitioning_allowed{true}, parent_mesh(&m){
50  if (me_is_multithreaded_now()) prohibit_partitioning();
51  mark_region_changed();
52  }
53 
54  mesh_region::mesh_region(const dal::bit_vector &bv)
55  : p(std::make_shared<impl>()), id_(size_type(-2)), type_(size_type(-1)),
56  partitioning_allowed{true}, parent_mesh(nullptr){
57  if (me_is_multithreaded_now()) prohibit_partitioning();
58  add(bv);
59  mark_region_changed();
60  }
61 
62  void mesh_region::mark_region_changed() const{
63  index_updated.all_threads() = false;
64  partitions_updated.all_threads() = false;
65  serial_index_updated = false;
66  }
67 
68  void mesh_region::touch_parent_mesh(){
69  if (parent_mesh) parent_mesh->touch_from_region(id_);
70  }
71 
72  const mesh_region& mesh_region::from_mesh(const mesh &m) const{
73  if (!p)
74  {
75  auto r = const_cast<mesh_region*>(this);
76  if (id_ == size_type(-1))
77  {
78  r->p = std::make_shared<impl>();
79  r->add(m.convex_index());
80  }
81  else if (id_ != size_type(-2))
82  {
83  *r = m.region(id_);
84  }
85  }
86  mark_region_changed();
87  return *this;
88  }
89 
90  mesh_region& mesh_region::operator=(const mesh_region &from){
91  if (!parent_mesh && !from.parent_mesh){
92  id_ = from.id_;
93  type_ = from.type_;
94  partitioning_allowed.store(from.partitioning_allowed.load());
95  if (from.p) {
96  if (!p) p = std::make_shared<impl>();
97  wp() = from.rp();
98  }
99  else p = nullptr;
100  }
101  else if (!parent_mesh) {
102  p = from.p;
103  id_ = from.id_;
104  type_ = from.type_;
105  parent_mesh = from.parent_mesh;
106  partitioning_allowed.store(from.partitioning_allowed.load());
107  }
108  else {
109  if (from.p){
110  wp() = from.rp();
111  type_= from.get_type();
112  partitioning_allowed.store(from.partitioning_allowed.load());
113  }
114  else if (from.id_ == size_type(-1)) {
115  clear();
116  add(parent_mesh->convex_index());
117  type_ = size_type(-1);
118  partitioning_allowed = true;
119  }
120  touch_parent_mesh();
121  }
122  mark_region_changed();
123  return *this;
124  }
125 
126  bool mesh_region::compare(const mesh &m1,
127  const mesh_region &mr,
128  const mesh &m2) const {
129  if (&m1 != &m2) return false;
130  if (!p && !mr.p) return (id_ == mr.id_);
131  from_mesh(m1);
132  mr.from_mesh(m2);
133  if (p && !(mr.p)) return false;
134  if (!p && mr.p) return false;
135  if (p)
136  if (p->m != mr.p->m) return false;
137  return true;
138  }
139 
140  face_bitset mesh_region::operator[](size_t cv) const{
141  auto it = rp().m.find(cv);
142  if (it != rp().m.end()) return it->second;
143  else return {};
144  }
145 
146  void mesh_region::update_partition_iterators() const{
147  if ((partitions_updated.thrd_cast() == true)) return;
148  itbegin = partition_begin();
149  itend = partition_end ();
150  partitions_updated = true;
151  }
152  mesh_region::const_iterator
153  mesh_region::partition_begin( ) const{
154  auto region_size = rp().m.size();
155  if (region_size < partitions_updated.num_threads()){
156  //for small regions: put the whole region into zero thread
157  if (partitions_updated.this_thread() == 0) return rp().m.begin();
158  else return rp().m.end();
159  }
160  auto partition_size = static_cast<size_type>
161  (std::ceil(static_cast<scalar_type>(region_size)/
162  static_cast<scalar_type >(partitions_updated.num_threads())));
163  auto index_begin = partition_size * partitions_updated.this_thread();
164  if (index_begin >= region_size ) return rp().m.end();
165 
166  auto it = rp().m.begin();
167  for (size_type i = 0; i != index_begin; ++i) ++it;
168  return it;
169  }
170 
171  mesh_region::const_iterator
172  mesh_region::partition_end( ) const{
173  auto region_size = rp().m.size();
174  if (region_size< partitions_updated.num_threads()) return rp().m.end();
175 
176  auto partition_size = static_cast<size_type>
177  (std::ceil(static_cast<scalar_type>(region_size)/
178  static_cast<scalar_type >(partitions_updated.num_threads())));
179  auto index_end = partition_size * (partitions_updated.this_thread() + 1);
180  if (index_end >= region_size ) return rp().m.end();
181 
182  auto it = rp().m.begin();
183  for (size_type i = 0; i < index_end && it != rp().m.end(); ++i) ++it;
184  return it;
185  }
186 
187  mesh_region::const_iterator mesh_region::begin() const{
188  GMM_ASSERT1(p != 0, "Internal error");
189  if (me_is_multithreaded_now() && partitioning_allowed){
190  update_partition_iterators();
191  return itbegin;
192  }
193  else return rp().m.begin();
194  }
195 
196  mesh_region::const_iterator mesh_region::end() const{
197  if (me_is_multithreaded_now() && partitioning_allowed){
198  update_partition_iterators();
199  return itend;
200  }
201  else return rp().m.end();
202  }
203 
205  partitioning_allowed = true;
206  }
207 
208  void mesh_region::bounding_box(base_node& Pmin, base_node& Pmax) const{
209  auto &mesh = *get_parent_mesh();
210  for (auto &&cv : dal::bv_iterable_c{index()}) {
211  for (const auto &pt : mesh.points_of_convex(cv)) {
212  for (auto j = 0; j != Pmin.size(); ++j){
213  Pmin[j] = std::min(Pmin[j], pt[j]);
214  Pmax[j] = std::max(Pmax[j], pt[j]);
215  }
216  }
217  }
218  }
219 
221  partitioning_allowed = false;
222  }
223 
224  bool mesh_region::is_partitioning_allowed() const{
225  return partitioning_allowed;
226  }
227 
228  void mesh_region::update_index() const{
229  auto& convex_index = me_is_multithreaded_now() ?
230  rp().index_.thrd_cast() : rp().serial_index_;
231  if (convex_index.card() != 0) convex_index.clear();
232  for (auto &&pair : *this){
233  if (pair.second.any()) convex_index.add(pair.first);
234  }
235  }
236 
237  const dal::bit_vector& mesh_region::index() const{
238  GMM_ASSERT1(p, "Use from_mesh on that region before");
239  if (me_is_multithreaded_now()) {
240  if (!(index_updated.thrd_cast() == true)){
241  update_index();
242  index_updated = true;
243  }
244  return rp().index_;
245  }
246  else {
247  if (!serial_index_updated){
248  update_index();
249  serial_index_updated = true;
250  }
251  return rp().serial_index_;
252  }
253  }
254 
255  void mesh_region::add(const dal::bit_vector &bv){
256  for (dal::bv_visitor i(bv); !i.finished(); ++i){
257  wp().m[i].set(0, 1);
258  }
259  touch_parent_mesh();
260  mark_region_changed();
261  }
262 
263  void mesh_region::add(size_type cv, short_type f){
264  wp().m[cv].set(short_type(f + 1), 1);
265  touch_parent_mesh();
266  mark_region_changed();
267  }
268 
269  void mesh_region::sup_all(size_type cv){
270  auto it = wp().m.find(cv);
271  if (it != wp().m.end()){
272  wp().m.erase(it);
273  touch_parent_mesh();
274  mark_region_changed();
275  }
276  }
277 
278  void mesh_region::sup(size_type cv, short_type f){
279  auto it = wp().m.find(cv);
280  if (it != wp().m.end()) {
281  it->second.set(short_type(f + 1), 0);
282  if (it->second.none()) wp().m.erase(it);
283  touch_parent_mesh();
284  mark_region_changed();
285  }
286  }
287 
288  void mesh_region::clear(){
289  wp().m.clear();
290  touch_parent_mesh();
291  mark_region_changed();
292  }
293 
294  void mesh_region::clean(){
295  for (map_t::iterator it = wp().m.begin(), itn;
296  it != wp().m.end(); it = itn) {
297  itn = it;
298  ++itn;
299  if (!(*it).second.any()) wp().m.erase(it);
300  }
301  touch_parent_mesh();
302  mark_region_changed();
303  }
304 
305  void mesh_region::swap_convex(size_type cv1, size_type cv2){
306  auto it1 = wp().m.find(cv1);
307  auto it2 = wp().m.find(cv2);
308  auto ite = wp().m.end();
309  face_bitset f1, f2;
310 
311  if (it1 != ite) f1 = it1->second;
312  if (it2 != ite) f2 = it2->second;
313  if (!f1.none()) wp().m[cv2] = f1;
314  else if (it2 != ite) wp().m.erase(it2);
315  if (!f2.none()) wp().m[cv1] = f2;
316  else if (it1 != ite) wp().m.erase(it1);
317  touch_parent_mesh();
318  mark_region_changed();
319  }
320 
321  bool mesh_region::is_in(size_type cv, short_type f) const{
322  GMM_ASSERT1(p, "Use from mesh on that region before");
323  auto it = rp().m.find(cv);
324  if (it == rp().m.end() || short_type(f+1) >= MAX_FACES_PER_CV) return false;
325  return ((*it).second)[short_type(f+1)];
326  }
327 
328  bool mesh_region::is_in(size_type cv, short_type f, const mesh &m) const{
329  if (p) {
330  auto it = rp().m.find(cv);
331  if (it == rp().m.end() || short_type(f+1) >= MAX_FACES_PER_CV)
332  return false;
333  return ((*it).second)[short_type(f+1)];
334  }
335  else{
336  if (id() == size_type(-1)) return true;
337  else return m.region(id()).is_in(cv, f);
338  }
339  }
340 
341  bool mesh_region::is_empty() const{
342  return rp().m.empty();
343  }
344 
346  return is_empty() ||
347  (or_mask()[0] == true && or_mask().count() == 1);
348  }
349 
351  return (id() != size_type(-1)) && (is_empty() || (and_mask()[0] == false));
352  }
353 
354  face_bitset mesh_region::faces_of_convex(size_type cv) const{
355  auto it = rp().m.find(cv);
356  if (it != rp().m.end()) return ((*it).second) >> 1;
357  else return face_bitset();
358  }
359 
360  face_bitset mesh_region::and_mask() const{
361  face_bitset bs;
362  if (rp().m.empty()) return bs;
363  bs.set();
364  for (auto it = rp().m.begin(); it != rp().m.end(); ++it)
365  if ( (*it).second.any() ) bs &= (*it).second;
366  return bs;
367  }
368 
369  face_bitset mesh_region::or_mask() const{
370  face_bitset bs;
371  if (rp().m.empty()) return bs;
372  for (auto it = rp().m.begin(); it != rp().m.end(); ++it)
373  if ( (*it).second.any() ) bs |= (*it).second;
374  return bs;
375  }
376 
378  size_type sz = 0;
379  for (auto it = begin(); it != end(); ++it)
380  sz += (*it).second.count();
381  return sz;
382  }
383 
384  size_type mesh_region::unpartitioned_size() const{
385  size_type sz = 0;
386  for (auto it = rp().m.begin(); it != rp().m.end(); ++it)
387  sz += (*it).second.count();
388  return sz;
389  }
390 
392  const mesh_region &b){
393  GMM_TRACE4("intersection of "<<a.id()<<" and "<<b.id());
394  mesh_region r;
395  /* we do not allow the "all_convexes" kind of regions
396  for these operations as there are not intended to be manipulated
397  (they only exist to provide a default argument to the mesh_region
398  parameters of assembly procedures etc. */
399  GMM_ASSERT1(a.id() != size_type(-1)||
400  b.id() != size_type(-1), "the 'all_convexes' regions "
401  "are not supported for set operations");
402  if (a.id() == size_type(-1)){
403  for (const_iterator it = b.begin();it != b.end(); ++it) r.wp().m.insert(*it);
404  return r;
405  }
406  else if (b.id() == size_type(-1)){
407  for (const_iterator it = a.begin();it != a.end(); ++it) r.wp().m.insert(*it);
408  return r;
409  }
410 
411  auto ita = a.begin(), enda = a.end(),
412  itb = b.begin(), endb = b.end();
413 
414  while (ita != enda && itb != endb) {
415  if (ita->first < itb->first) ++ita;
416  else if (ita->first > itb->first) ++itb;
417  else {
418  face_bitset maska = ita->second, maskb = itb->second, bs;
419  if (maska[0] && !maskb[0]) bs = maskb;
420  else if (maskb[0] && !maska[0]) bs = maska;
421  else bs = maska & maskb;
422  if (bs.any()) r.wp().m.insert(r.wp().m.end(), std::make_pair(ita->first,bs));
423  ++ita; ++itb;
424  }
425  }
426  return r;
427  }
428 
430  const mesh_region &b){
431  GMM_TRACE4("Merger of " << a.id() << " and " << b.id());
432  mesh_region r;
433  GMM_ASSERT1(a.id() != size_type(-1) &&
434  b.id() != size_type(-1), "the 'all_convexes' regions "
435  "are not supported for set operations");
436  for (auto it = a.begin();it != a.end(); ++it){
437  r.wp().m.insert(*it);
438  }
439  for (auto it = b.begin();it != b.end(); ++it){
440  r.wp().m[it->first] |= it->second;
441  }
442  return r;
443  }
444 
445 
447  const mesh_region &b){
448  GMM_TRACE4("subtraction of "<<a.id()<<" and "<<b.id());
449  mesh_region r;
450  GMM_ASSERT1(a.id() != size_type(-1) &&
451  b.id() != size_type(-1), "the 'all_convexes' regions "
452  "are not supported for set operations");
453  for (auto ita = a.begin();ita != a.end(); ++ita)
454  r.wp().m.insert(*ita);
455 
456  for (auto itb = b.begin();itb != b.end(); ++itb){
457  auto cv = itb->first;
458  auto it = r.wp().m.find(cv);
459  if (it != r.wp().m.end()){
460  it->second &= ~(itb->second);
461  if (it->second.none()) r.wp().m.erase(it);
462  }
463  }
464  return r;
465  }
466 
468  const mesh_region &rg2,
469  const getfem::mesh& m2) const{
470  if (&m1 != &m2) return 0;
471  int r = 1, partially = 0;
472  for (mr_visitor cv(*this, m1); !cv.finished(); cv.next())
473  if (cv.is_face() && rg2.is_in(cv.cv(),short_type(-1), m2))
474  partially = -1;
475  else
476  r = 0;
477  return r == 1 ? 1 : partially;
478  }
479 
481  return m.regions_index().last_true()+1;
482  }
483 
484  void mesh_region::error_if_not_faces() const{
485  GMM_ASSERT1(is_only_faces(), "Expecting a set of faces, not convexes");
486  }
487 
488  void mesh_region::error_if_not_convexes() const{
489  GMM_ASSERT1(is_only_convexes(), "Expecting a set of convexes, not faces");
490  }
491 
492  void mesh_region::error_if_not_homogeneous() const{
493  GMM_ASSERT1(is_only_faces() || is_only_convexes(), "Expecting a set "
494  "of convexes or a set of faces, but not a mixed set");
495  }
496 
497 
498 #if GETFEM_PARA_LEVEL > 1
499 
500  mesh_region::visitor::visitor(const mesh_region &s, const mesh &m,
501  bool intersect_with_mpi) :
502  cv_(size_type(-1)), f_(short_type(-1)), finished_(false)
503  {
504  if ((me_is_multithreaded_now() && s.partitioning_allowed)) {
505  s.from_mesh(m);
506  init(s);
507  } else {
508  if (s.id() == size_type(-1)) {
509  if (intersect_with_mpi)
510  init(m.get_mpi_region());
511  else
512  init(m.convex_index());
513  } else if (s.id() == size_type(-2) && s.p.get()) {
514  if (intersect_with_mpi) {
515  mpi_rg = std::make_unique<mesh_region>(s);
516  mpi_rg->from_mesh(m);
517  m.intersect_with_mpi_region(*mpi_rg);
518  init(*mpi_rg);
519  } else
520  init(s);
521  } else {
522  GMM_ASSERT1(s.id() != size_type(-2), "Internal error");
523  if (intersect_with_mpi)
524  init(m.get_mpi_sub_region(s.id()));
525  else
526  init(m.region(s.id()));
527  }
528  }
529  }
530 
531 #else
532 
533  mesh_region::visitor::visitor(const mesh_region &s, const mesh &m, bool)
534  :cv_(size_type(-1)), f_(short_type(-1)), finished_(false){
535  if ((me_is_multithreaded_now() && s.partitioning_allowed)) {
536  s.from_mesh(m);
537  init(s);
538  }
539  else {
540  if (s.id() == size_type(-1)) {
541  init(m.convex_index());
542  } else if (s.id() == size_type(-2) && s.p.get()) {
543  init(s);
544  } else {
545  GMM_ASSERT1(s.id() != size_type(-2), "Internal error");
546  init(m.region(s.id()));
547  }
548  }
549  }
550 
551 #endif
552 
553  bool mesh_region::visitor::next(){
554  if (whole_mesh) {
555  if (itb == iteb) {
556  finished_ = true;
557  return false;
558  }
559  cv_ = itb.index();
560  c = 0;
561  f_ = 0;
562  ++itb;
563  while (itb != iteb && !(*itb)) ++itb;
564  return true;
565  }
566  while (c.none()){
567  if (it == ite) { finished_=true; return false; }
568  cv_ = it->first;
569  c = it->second;
570  f_ = short_type(-1);
571  ++it;
572  if (c.none()) continue;
573  }
574  next_face();
575  return true;
576  }
577 
578  mesh_region::visitor::visitor(const mesh_region &s) :
579  cv_(size_type(-1)), f_(short_type(-1)), finished_(false){
580  init(s);
581  }
582 
583  void mesh_region::visitor::init(const dal::bit_vector &bv){
584  whole_mesh = true;
585  itb = bv.begin(); iteb = bv.end();
586  while (itb != iteb && !(*itb)) ++itb;
587  next();
588  }
589 
590  void mesh_region::visitor::init(const mesh_region &s){
591  whole_mesh = false;
592  it = s.begin();
593  ite = s.end();
594  next();
595  }
596 
597  std::ostream & operator <<(std::ostream &os, const mesh_region &w){
598  if (w.id() == size_type(-1))
599  os << " ALL_CONVEXES";
600  else if (w.p){
601  for (mr_visitor cv(w); !cv.finished(); cv.next()){
602  os << cv.cv();
603  if (cv.is_face()) os << "/" << cv.f();
604  os << " ";
605  }
606  }
607  else{
608  os << " region " << w.id();
609  }
610  return os;
611  }
612 
613  struct dummy_mesh_region_{
614  mesh_region mr;
615  };
616 
619  }
620 
621 } /* end of namespace getfem. */
dal::singleton::instance
static T & instance()
Instance from the current thread.
Definition: dal_singleton.h:165
getfem::me_is_multithreaded_now
bool me_is_multithreaded_now()
is the program running in the parallel section
Definition: getfem_omp.cc:106
getfem_mesh_region.h
region objects (set of convexes and/or convex faces)
getfem::mesh::region
const mesh_region region(size_type id) const
Return the region of index 'id'.
Definition: getfem_mesh.h:421
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
getfem::mesh_region::subtract
static mesh_region subtract(const mesh_region &a, const mesh_region &b)
subtract the second region from the first one
Definition: getfem_mesh_region.cc:446
getfem::mesh_region::size
size_type size() const
Region size, or the size of the region partition on the current thread if the region is partitioned.
Definition: getfem_mesh_region.cc:377
getfem::mesh_region::region_is_faces_of
int region_is_faces_of(const getfem::mesh &m1, const mesh_region &rg2, const getfem::mesh &m2) const
Test if the region is a boundary of a list of faces of elements of region rg.
Definition: getfem_mesh_region.cc:467
getfem::dummy_mesh_region
const mesh_region & dummy_mesh_region()
Dummy mesh_region for default parameter of functions.
Definition: getfem_mesh_region.cc:617
getfem::mesh_region::allow_partitioning
void allow_partitioning()
In multithreaded part of the program makes only a partition of the region visible in the index() and ...
Definition: getfem_mesh_region.cc:204
getfem::mesh_region::prohibit_partitioning
void prohibit_partitioning()
Disregard partitioning, which means being able to see the whole region in multithreaded code.
Definition: getfem_mesh_region.cc:220
getfem_mesh.h
Define a getfem::getfem_mesh object.
getfem::mesh_region::free_region_id
static size_type free_region_id(const getfem::mesh &m)
Extract the next region number that does not yet exists in the mesh.
Definition: getfem_mesh_region.cc:480
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
getfem::mesh_region::bounding_box
void bounding_box(base_node &Pmin, base_node &Pmax) const
Return the bounding box [Pmin - Pmax] of the mesh_region.
Definition: getfem_mesh_region.cc:208
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_region::is_only_faces
bool is_only_faces() const
Return true if the region do contain only convex faces.
Definition: getfem_mesh_region.cc:350
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
getfem_omp.h
Tools for multithreaded, OpenMP and Boost based parallelization.
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_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
getfem::mesh
Describe a mesh (collection of convexes (elements) and points).
Definition: getfem_mesh.h:95
getfem::mesh_region
structure used to hold a set of convexes and/or convex faces.
Definition: getfem_mesh_region.h:55
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::mesh_region::merge
static mesh_region merge(const mesh_region &a, const mesh_region &b)
return the union of two mesh_regions
Definition: getfem_mesh_region.cc:429