GetFEM  5.4.2
getfem_mesh_slice.cc
1 /*===========================================================================
2 
3  Copyright (C) 2003-2020 Julien Pommier
4 
5  This file is a part of GetFEM
6 
7  GetFEM is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program; if not, write to the Free Software Foundation,
18  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
19 
20 ===========================================================================*/
21 
24 
25 namespace getfem {
26 
27  std::ostream& operator<<(std::ostream& o, const stored_mesh_slice& m) {
28  o << "stored_mesh_slice, containing " << m.nb_convex() << " convexes\n";
29  for (size_type ic = 0; ic < m.nb_convex(); ++ic) {
30  o << "slice convex #" << ic << " (original = " << m.convex_num(ic)
31  << ")\n";
32  for (size_type i = 0; i < m.nodes(ic).size(); ++i) {
33  o << "node " << i << ": " << m.nodes(ic)[i].pt << ", ref="
34  << m.nodes(ic)[i].pt_ref << " flist=" << m.nodes(ic)[i].faces
35  << endl;
36  }
37  for (size_type i = 0; i < m.simplexes(ic).size(); ++i) {
38  o << "simplex " << i << ", inodes=";
39  for (size_type j=0;j< m.simplexes(ic)[i].dim()+1;++j)
40  o << m.simplexes(ic)[i].inodes[j] << " ";
41  o << endl;
42  }
43  o << endl;
44  }
45  return o;
46  }
47 
48  void stored_mesh_slice::write_to_file(const std::string &name,
49  bool with_mesh) const {
50  std::ofstream o(name.c_str());
51  GMM_ASSERT1(o, "impossible to open file '" << name << "'");
52  o << "% GETFEM SLICE FILE " << '\n';
53  o << "% GETFEM VERSION " << GETFEM_VERSION << '\n' << '\n' << '\n';
54  if (with_mesh) linked_mesh().write_to_file(o);
55  write_to_file(o);
56  }
57 
58  void stored_mesh_slice::write_to_file(std::ostream &os) const {
59  os << "\nBEGIN MESH_SLICE\n";
60  os << " DIM " << int(dim()) << "\n";
61  for (unsigned i=0; i < cvlst.size(); ++i) {
62  const convex_slice &cs = cvlst[i];
63  os << " CONVEX " << cs.cv_num
64  << " " << int(cs.fcnt)
65  << " " << int(cs.discont) << "\n"
66  << " " << cs.nodes.size() << " " << cs.simplexes.size() << "\n";
67  for (unsigned j=0; j < cs.nodes.size(); ++j) {
68  os << "\t";
69  for (unsigned k=0; k < cs.nodes[j].pt.size(); ++k) {
70  if (k) os << " ";
71  os << cs.nodes[j].pt[k];
72  }
73  os << ";";
74  for (unsigned k=0; k < cs.nodes[j].pt_ref.size(); ++k)
75  os << " " << cs.nodes[j].pt_ref[k];
76  os << "; "; os << cs.nodes[j].faces.to_ulong();;
77  os << "\n";
78  }
79  for (unsigned j=0; j < cs.simplexes.size(); ++j) {
80  os << "\t" << cs.simplexes[j].inodes.size() << ":";
81  for (unsigned k=0; k < cs.simplexes[j].inodes.size(); ++k) {
82  os << " " << cs.simplexes[j].inodes[k];
83  }
84  os << "\n";
85  }
86  }
87  os << "END MESH_SLICE\n";
88  }
89 
90  void stored_mesh_slice::read_from_file(const std::string &name,
91  const getfem::mesh &m) {
92  std::ifstream o(name.c_str());
93  GMM_ASSERT1(o, "slice file '" << name << "' does not exist");
94  read_from_file(o,m);
95  }
96 
97  void stored_mesh_slice::read_from_file(std::istream &ist,
98  const getfem::mesh &m) {
99  if (!poriginal_mesh) {
100  poriginal_mesh = &m;
101  } else GMM_ASSERT1(poriginal_mesh == &m, "wrong mesh..");
102 
103  dim_ = m.dim();
104  cv2pos.clear();
105  cv2pos.resize(m.nb_allocated_convex(), size_type(-1));
106 
107  std::string tmp;
108  ist.precision(16);
109  ist.seekg(0);ist.clear();
110  bgeot::read_until(ist, "BEGIN MESH_SLICE");
111 
112  mesh_slicer::cs_nodes_ct nod;
113  mesh_slicer::cs_simplexes_ct sim;
114 
115 
116  while (true) {
117  ist >> std::ws; bgeot::get_token(ist, tmp);
118  if (bgeot::casecmp(tmp, "END")==0) {
119  break;
120  } else if (bgeot::casecmp(tmp, "DIM")==0) {
121  int d; ist >> d;
122  dim_ = d;
123  } else if (bgeot::casecmp(tmp, "CONVEX")==0) {
124  bgeot::get_token(ist,tmp);
125  size_type ic = atoi(tmp.c_str());
126  GMM_ASSERT1(m.convex_index().is_in(ic), "Convex " << ic <<
127  " does not exist, are you sure "
128  "that the mesh attached to this object is right one ?");
129  bgeot::pconvex_ref cvr = m.trans_of_convex(ic)->convex_ref();
130  unsigned fcnt, discont, nbn, nbs;
131  ist >> fcnt >> discont >> nbn >> nbs;
132  nod.resize(nbn);
133  sim.resize(nbs);
134  for (unsigned i=0; i < nbn; ++i) {
135  nod[i].pt.resize(dim());
136  nod[i].pt_ref.resize(cvr->structure()->dim());
137  for (unsigned j=0; j < dim(); ++j)
138  ist >> nod[i].pt[j];
139  ist >> bgeot::skip(";");
140  for (unsigned j=0; j < cvr->structure()->dim(); ++j)
141  ist >> nod[i].pt_ref[j];
142  ist >> bgeot::skip(";");
143  unsigned long ul; ist >> ul;
144  nod[i].faces = slice_node::faces_ct(int(ul));
145  }
146  for (unsigned i=0; i < nbs; ++i) {
147  unsigned np(0);
148  ist >> np >> bgeot::skip(":");
149  GMM_ASSERT1(np <= dim()+1, "invalid simplex..");
150  sim[i].inodes.resize(np);
151  for (unsigned j=0; j < np; ++j)
152  ist >> sim[i].inodes[j];
153  }
154  dal::bit_vector bv; bv.add(0, nbs);
155  set_convex(ic, cvr, nod, sim, dim_type(fcnt), bv, discont);
156  } else if (tmp.size()) {
157  GMM_ASSERT1(false, "Unexpected token '" << tmp <<
158  "' [pos=" << std::streamoff(ist.tellg()) << "]");
159  } else if (ist.eof()) {
160  GMM_ASSERT1(false, "Unexpected end of stream "
161  << "(missing BEGIN MESH_SLICE/END MESH_SLICE ?)");
162  }
163  }
164  }
165 
166  void slicer_build_stored_mesh_slice::exec(mesh_slicer &ms) {
167  if (!sl.poriginal_mesh) {
168  sl.poriginal_mesh = &ms.m;
169  sl.dim_ = sl.linked_mesh().dim();
170  sl.cv2pos.resize(sl.linked_mesh().nb_allocated_convex());
171  gmm::fill(sl.cv2pos, size_type(-1));
172  } else if (sl.poriginal_mesh != &ms.m) GMM_ASSERT1(false, "wrong mesh..");
173  sl.set_convex(ms.cv, ms.cvr, ms.nodes, ms.simplexes, dim_type(ms.fcnt),
174  ms.splx_in, ms.discont);
175  }
176 
177  void stored_mesh_slice::set_convex(size_type cv, bgeot::pconvex_ref cvr,
178  mesh_slicer::cs_nodes_ct cv_nodes,
179  mesh_slicer::cs_simplexes_ct cv_simplexes,
180  dim_type fcnt,
181  const dal::bit_vector& splx_in,
182  bool discont) {
183  /* push the used nodes and simplexes in the final list */
184  if (splx_in.card() == 0) return;
185  merged_nodes_available = false;
186  std::vector<size_type> nused(cv_nodes.size(), size_type(-1));
187  convex_slice *sc = 0;
188  GMM_ASSERT1(cv < cv2pos.size(), "internal error");
189  if (cv2pos[cv] == size_type(-1)) {
190  cv2pos[cv] = cvlst.size();
191  cvlst.push_back(convex_slice());
192  sc = &cvlst.back();
193  sc->cv_num = cv;
194  sc->cv_dim = cvr->structure()->dim();
195  sc->cv_nbfaces = dim_type(cvr->structure()->nb_faces());
196  sc->fcnt = fcnt;
197  sc->global_points_count = points_cnt;
198  sc->discont = discont;
199  } else {
200  sc = &cvlst[cv2pos[cv]];
201  assert(sc->cv_num == cv);
202  }
203  for (dal::bv_visitor snum(splx_in); !snum.finished(); ++snum) {
204  slice_simplex& s = cv_simplexes[snum];
205  for (size_type i=0; i < s.dim()+1; ++i) {
206  size_type lnum = s.inodes[i];
207  if (nused[lnum] == size_type(-1)) {
208  nused[lnum] = sc->nodes.size(); sc->nodes.push_back(cv_nodes[lnum]);
209  dim_ = std::max(int(dim_), int(cv_nodes[lnum].pt.size()));
210  points_cnt++;
211  }
212  s.inodes[i] = nused[lnum];
213  }
214  simplex_cnt.resize(dim_+1, 0);
215  simplex_cnt[cv_simplexes[snum].dim()]++;
216  sc->simplexes.push_back(cv_simplexes[snum]);
217  }
218  }
219 
220  struct get_edges_aux {
221  size_type iA, iB;
222  mutable bool slice_edge;
223  get_edges_aux(size_type a, size_type b, bool slice_edge_) :
224  iA(std::min(a,b)), iB(std::max(a,b)), slice_edge(slice_edge_) {}
225  bool operator<(const get_edges_aux& other) const {
226  /* ignore the slice_edge on purpose */
227  return (iA < other.iA || (iA == other.iA && iB < other.iB));
228  }
229  };
230 
231  void stored_mesh_slice::get_edges(std::vector<size_type> &edges,
232  dal::bit_vector &slice_edges,
233  bool from_merged_nodes) const {
234  if (from_merged_nodes && !merged_nodes_available) merge_nodes();
235  std::set<get_edges_aux> e;
236  for (cvlst_ct::const_iterator it=cvlst.begin(); it != cvlst.end(); ++it) {
237  for (size_type is=0; is < it->simplexes.size(); ++is) {
238  const slice_simplex &s = it->simplexes[is];
239  for (size_type i=0; i < s.dim(); ++i) {
240  for (size_type j=i+1; j <= s.dim(); ++j) {
241  const slice_node& A = it->nodes[s.inodes[i]];
242  const slice_node& B = it->nodes[s.inodes[j]];
243  /* duplicate with slicer_build_edges_mesh which also
244  builds a list of edges */
245  if ((A.faces & B.faces).count() >= unsigned(it->cv_dim-1)) {
246  slice_node::faces_ct fmask((1 << it->cv_nbfaces)-1);
247  fmask.flip();
248  size_type iA, iB;
249  iA = it->global_points_count + s.inodes[i];
250  iB = it->global_points_count + s.inodes[j];
251  if (from_merged_nodes) {
252  iA = to_merged_index[iA]; iB = to_merged_index[iB];
253  }
254  get_edges_aux a(iA,iB,((A.faces & B.faces) & fmask).any());
255  std::set<get_edges_aux>::iterator p=e.find(a);
256  if (p != e.end()) {
257  if (p->slice_edge && !a.slice_edge) p->slice_edge = false;
258  } else e.insert(a);
259  }
260  }
261  }
262  }
263  }
264  slice_edges.clear(); slice_edges.sup(0, e.size());
265  edges.clear(); edges.reserve(2*e.size());
266  for (std::set<get_edges_aux>::const_iterator p=e.begin();
267  p != e.end(); ++p) {
268  if (p->slice_edge) slice_edges.add(edges.size()/2);
269  edges.push_back(p->iA);edges.push_back(p->iB);
270  }
271  }
272 
273  void stored_mesh_slice::build(const getfem::mesh& m,
274  const slicer_action *a,
275  const slicer_action *b,
276  const slicer_action *c,
277  size_type nrefine) {
278  clear();
279  mesh_slicer slicer(m);
280  slicer.push_back_action(*const_cast<slicer_action*>(a));
281  if (b) slicer.push_back_action(*const_cast<slicer_action*>(b));
282  if (c) slicer.push_back_action(*const_cast<slicer_action*>(c));
283  slicer_build_stored_mesh_slice sbuild(*this);
284  slicer.push_back_action(sbuild);
285  slicer.exec(nrefine);
286  }
287 
288  void stored_mesh_slice::replay(slicer_action *a, slicer_action *b,
289  slicer_action *c) const {
290  mesh_slicer slicer(linked_mesh());
291  slicer.push_back_action(*a);
292  if (b) slicer.push_back_action(*b);
293  if (c) slicer.push_back_action(*c);
294  slicer.exec(*this);
295  }
296 
298  dim_ = newdim;
299  for (size_type ic=0; ic < nb_convex(); ++ic) {
300  for (mesh_slicer::cs_nodes_ct::iterator it=nodes(ic).begin();
301  it != nodes(ic).end(); ++it) {
302  it->pt.resize(newdim);
303  }
304  }
305  }
306 
308  GMM_ASSERT1(dim()==sl.dim(), "inconsistent dimensions for slice merging");
309  clear_merged_nodes();
310  cv2pos.resize(std::max(cv2pos.size(), sl.cv2pos.size()), size_type(-1));
311  for (size_type i=0; i < sl.nb_convex(); ++i)
312  GMM_ASSERT1(cv2pos[sl.convex_num(i)] == size_type(-1) ||
313  cvlst[cv2pos[sl.convex_num(i)]].cv_dim == sl.cvlst[i].cv_num,
314  "inconsistent dimensions for convex " << sl.cvlst[i].cv_num
315  << " on the slices");
316  for (size_type i=0; i < sl.nb_convex(); ++i) {
317  size_type cv = sl.convex_num(i);
318  if (cv2pos[cv] == size_type(-1)) {
319  cv2pos[cv] = cvlst.size();
320  cvlst.push_back(convex_slice());
321  }
322  const stored_mesh_slice::convex_slice *src = &sl.cvlst[i];
323  stored_mesh_slice::convex_slice *dst = &cvlst[cv2pos[cv]];
324  size_type n = dst->nodes.size();
325  dst->nodes.insert(dst->nodes.end(), src->nodes.begin(),
326  src->nodes.end());
327  for (mesh_slicer::cs_simplexes_ct::const_iterator
328  it = src->simplexes.begin(); it != src->simplexes.end(); ++it) {
329  dst->simplexes.push_back(*it);
330  for (size_type j = 0; j < (*it).dim()+1; ++j)
331  dst->simplexes.back().inodes[j] += n;
332  simplex_cnt[dst->simplexes.back().dim()]++;
333  }
334  points_cnt += src->nodes.size();
335  }
336  size_type count = 0;
337  for (size_type ic=0; ic < nb_convex(); ++ic) {
338  cvlst[ic].global_points_count = count; count += nodes(ic).size();
339  }
340  assert(count == points_cnt);
341  }
342 
343  void stored_mesh_slice::clear_merged_nodes() const {
344  merged_nodes_idx.clear(); merged_nodes.clear();
345  to_merged_index.clear();
346  merged_nodes_available = false;
347  }
348 
350  size_type count = 0;
351  mesh mp;
352  clear_merged_nodes();
353  std::vector<size_type> iv;
354  std::vector<const slice_node*> nv(nb_points());
355  to_merged_index.resize(nb_points());
356  for (cvlst_ct::const_iterator it = cvlst.begin(); it != cvlst.end(); ++it) {
357  for (size_type i=0; i < it->nodes.size(); ++i) {
358  nv[count] = &it->nodes[i];
359  to_merged_index[count++] = mp.add_point(it->nodes[i].pt);
360  // cout << "orig[" << count-1 << "] = " << nv[count-1]->pt
361  // << ", idx=" << to_merged_index[count-1] << "\n";
362  }
363  }
364  gmm::sorted_indexes(to_merged_index,iv);
365  //cout << "to_merged_index = " << iv << "\n";
366  merged_nodes.resize(nb_points());
367  merged_nodes_idx.reserve(nb_points()/8);
368 
369  merged_nodes_idx.push_back(0);
370  for (size_type i=0; i < nb_points(); ++i) {
371  merged_nodes[i].P = nv[iv[i]];
372  merged_nodes[i].pos = unsigned(iv[i]);
373  // cout << "i=" << i << " -> {" << merged_nodes[i].P->pt
374  // << "," << merged_nodes[i].pos << "}\n";
375  if (i == nb_points()-1 ||
376  to_merged_index[iv[i+1]] != to_merged_index[iv[i]])
377  merged_nodes_idx.push_back(i+1);
378  }
379  //cout << "merged_nodes_idx = " << merged_nodes_idx << "\n";
380  merged_nodes_available = true;
381  }
382 
383  size_type stored_mesh_slice::memsize() const {
384  size_type sz = sizeof(stored_mesh_slice);
385  for (cvlst_ct::const_iterator it = cvlst.begin();
386  it != cvlst.end(); ++it) {
387  sz += sizeof(size_type);
388  // cerr << "memsize: convex " << it->cv_num << " nodes:"
389  // << it->nodes.size() << ", splxs:" << it->simplexes.size()
390  // << ", sz=" << sz << "\n";
391  for (size_type i=0; i < it->nodes.size(); ++i) {
392  // cerr << " point " << i << ": size+= " << sizeof(slice_node)
393  // << "+" << it->nodes[i].pt.memsize() << "+"
394  // << it->nodes[i].pt_ref.memsize() << "-"
395  // << sizeof(it->nodes[i].pt)*2 << "\n";*/
396  sz += sizeof(slice_node) +
397  (it->nodes[i].pt.memsize()+it->nodes[i].pt_ref.memsize())
398  - sizeof(it->nodes[i].pt)*2;
399  }
400  for (size_type i=0; i < it->simplexes.size(); ++i) {
401  // cerr << " simplex " << i << ": size+= " << sizeof(slice_simplex)
402  // << "+" << it->simplexes[i].inodes.size()*sizeof(size_type)
403  // << "\n";*/
404  sz += sizeof(slice_simplex) +
405  it->simplexes[i].inodes.size()*sizeof(size_type);
406  }
407  }
408  sz += cv2pos.size() * sizeof(size_type);
409  return sz;
410  }
411 
412 }
getfem::slicer_action
generic slicer class.
Definition: getfem_mesh_slicers.h:237
bgeot::operator<<
std::ostream & operator<<(std::ostream &o, const convex_structure &cv)
Print the details of the convex structure cvs to the output stream o.
Definition: bgeot_convex_structure.cc:71
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
getfem::slicer_build_stored_mesh_slice
a getfem::mesh_slicer whose side effect is to build a stored_mesh_slice object.
Definition: getfem_mesh_slice.h:303
getfem::stored_mesh_slice
The output of a getfem::mesh_slicer which has been recorded.
Definition: getfem_mesh_slice.h:47
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
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
getfem::mesh_slicer
Apply a serie a slicing operations to a mesh.
Definition: getfem_mesh_slicers.h:103
getfem::stored_mesh_slice::write_to_file
void write_to_file(std::ostream &os) const
Save a slice content to a text file.
Definition: getfem_mesh_slice.cc:58
getfem::stored_mesh_slice::get_edges
void get_edges(std::vector< size_type > &edges, dal::bit_vector &slice_edges, bool from_merged_nodes) const
Extract the list of mesh edges.
Definition: getfem_mesh_slice.cc:231
getfem::stored_mesh_slice::nb_points
size_type nb_points() const
Return the number of nodes in the slice.
Definition: getfem_mesh_slice.h:118
getfem::stored_mesh_slice::set_dim
void set_dim(size_type newdim)
change the slice dimension (append zeros or truncate node coordinates..)
Definition: getfem_mesh_slice.cc:297
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::stored_mesh_slice::merge
void merge(const stored_mesh_slice &sl)
merge with another mesh slice.
Definition: getfem_mesh_slice.cc:307
getfem::stored_mesh_slice::nb_convex
size_type nb_convex() const
return the number of convexes of the original mesh referenced in the slice
Definition: getfem_mesh_slice.h:98
getfem::stored_mesh_slice::merge_nodes
void merge_nodes() const
build a list of merged nodes.
Definition: getfem_mesh_slice.cc:349
getfem_mesh_slice.h
Define the class getfem::stored_mesh_slice.
getfem::stored_mesh_slice::read_from_file
void read_from_file(std::istream &ist, const getfem::mesh &m)
Read a slice from a file.
Definition: getfem_mesh_slice.cc:97
bgeot_geotrans_inv.h
Inversion of geometric transformations.
getfem::stored_mesh_slice::dim
size_type dim() const
return the slice dimension
Definition: getfem_mesh_slice.h:108
getfem::stored_mesh_slice::nodes
const mesh_slicer::cs_nodes_ct & nodes(size_type ic) const
Return the list of nodes for the 'ic'th convex of the slice.
Definition: getfem_mesh_slice.h:120
getfem::stored_mesh_slice::build
void build(const getfem::mesh &m, const slicer_action &a, size_type nrefine=1)
Build the slice, by applying a slicer_action operation.
Definition: getfem_mesh_slice.h:197
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::mesh
Describe a mesh (collection of convexes (elements) and points).
Definition: getfem_mesh.h:95
getfem::stored_mesh_slice::convex_num
size_type convex_num(size_type ic) const
return the original convex number of the 'ic'th convex referenced in the slice
Definition: getfem_mesh_slice.h:101
getfem::stored_mesh_slice::replay
void replay(slicer_action &a) const
Apply the listed slicer_action(s) to the slice object.
Definition: getfem_mesh_slice.h:215
getfem::stored_mesh_slice::linked_mesh
const mesh & linked_mesh() const
return a pointer to the original mesh
Definition: getfem_mesh_slice.h:110