v0.14.0
TetGenInterface.cpp
Go to the documentation of this file.
1 /** \file TetGenInterface.cpp
2  * \brief TetGen interface for meshing/re-meshing and on the fly mesh creation
3  *
4  */
5 
6 
7 
8 #ifdef WITH_TETGEN
9 #include <tetgen.h>
10 #endif
11 
12 #ifdef WITH_TETGEN
13 
14 // #define DEBUG_TETGEN
15 #ifdef DEBUG_TETGEN
16 #include <FTensor.hpp>
17 
18 template <class T1, class T2>
19 static inline MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det) {
21  det = +t(0, 0) * t(1, 1) * t(2, 2) + t(1, 0) * t(2, 1) * t(0, 2) +
22  t(2, 0) * t(0, 1) * t(1, 2) - t(0, 0) * t(2, 1) * t(1, 2) -
23  t(2, 0) * t(1, 1) * t(0, 2) - t(1, 0) * t(0, 1) * t(2, 2);
25 }
26 #endif
27 
28 namespace MoFEM {
29 
31 TetGenInterface::query_interface(boost::typeindex::type_index type_index,
32  UnknownInterface **iface) const {
33  *iface = const_cast<TetGenInterface *>(this);
34  return 0;
35 }
36 
38 TetGenInterface::inData(Range &ents, tetgenio &in,
39  std::map<EntityHandle, unsigned long> &moab_tetgen_map,
40  std::map<unsigned long, EntityHandle> &tetgen_moab_map,
41  Tag th) {
43 
44  Interface &m_field = cOre;
45  Range::iterator it;
46 
47  Tag th_marker;
48  int def_marker = 0;
49  CHKERR m_field.get_moab().tag_get_handle(
50  "TETGEN_MARKER", 1, MB_TYPE_INTEGER, th_marker,
51  MB_TAG_CREAT | MB_TAG_SPARSE, &def_marker);
52 
53  // All indices start from 0
54  in.firstnumber = 0;
55 
56  Range points = ents.subset_by_dimension(0);
57  in.numberofpoints = points.size();
58  if (points.size() > 0) {
59  in.pointlist = new double[in.numberofpoints * 3];
60  in.pointmarkerlist = new int[in.numberofpoints];
61  if (th) {
62  CHKERR m_field.get_moab().tag_get_data(th, points, in.pointlist);
63  } else {
64  CHKERR m_field.get_moab().get_coords(points, in.pointlist);
65  }
66  CHKERR m_field.get_moab().tag_get_data(th_marker, points,
67  in.pointmarkerlist);
68  it = points.begin();
69  for (int ii = 0; it != points.end(); it++, ii++) {
70  unsigned long iii = MBVERTEX | (ii << MB_TYPE_WIDTH);
71  tetgen_moab_map[iii] = *it;
72  moab_tetgen_map[*it] = iii;
73  }
74  }
75 
76  Range tets = ents.subset_by_type(MBTET);
77  in.numberoftetrahedra = tets.size();
78  if (in.numberoftetrahedra > 0) {
79  in.tetrahedronlist = new int[4 * ents.subset_by_type(MBTET).size()];
80  it = tets.begin();
81  for (int ii = 0; it != tets.end(); it++, ii++) {
82  int num_nodes;
83  const EntityHandle *conn;
84  CHKERR m_field.get_moab().get_connectivity(*it, conn, num_nodes, true);
85 #ifdef DEBUG_TETGEN
86  {
87  double coords[12];
88  if (th) {
89  CHKERR m_field.tag_get_data(th, conn, num_nodes, coords);
90  } else {
91  CHKERR m_field.get_moab().get_coords(conn, num_nodes, coords);
92  }
93  double diff_n[12];
94  ShapeDiffMBTET(diff_n);
95  FTensor::Tensor1<double *, 3> t_diff_n(&diff_n[0], &diff_n[1],
96  &diff_n[2], 3);
97  FTensor::Tensor1<double *, 3> t_coords(&coords[0], &coords[1],
98  &coords[2], 3);
102  jac(i, j) = 0;
103  for (int nn = 0; nn != 4; nn++) {
104  jac(i, j) += t_coords(i) * t_diff_n(j);
105  ++t_coords;
106  ++t_diff_n;
107  }
108  double det;
109  determinantTensor3by3(jac, det);
110  if (det <= 0) {
111  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
112  "Negative volume det = %6.4e", det);
113  }
114  }
115 #endif
116  tetgen_moab_map[MBTET | (ii << MB_TYPE_WIDTH)] = *it;
117  moab_tetgen_map[*it] = MBTET | (ii << MB_TYPE_WIDTH);
118  for (int nn = 0; nn != 4; nn++) {
119  if (moab_tetgen_map.find(conn[nn]) == moab_tetgen_map.end()) {
120  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
121  "data inconsistency between TetGen and MoAB");
122  }
123  in.tetrahedronlist[4 * ii + nn] =
124  moab_tetgen_map[conn[nn]] >> MB_TYPE_WIDTH;
125  }
126  }
127  }
128 
129  Range tris = ents.subset_by_type(MBTRI);
130  in.numberoftrifaces = tris.size();
131  if (in.numberoftrifaces) {
132  in.trifacelist = new int[3 * in.numberoftrifaces];
133  in.trifacemarkerlist = new int[in.numberoftrifaces];
134  // std::fill(&in.trifacemarkerlist[0],&in.trifacemarkerlist[in.numberoftrifaces],1);
135  CHKERR m_field.get_moab().tag_get_data(th_marker, tris,
136  in.trifacemarkerlist);
137  it = tris.begin();
138  for (int ii = 0; it != tris.end(); it++, ii++) {
139  int num_nodes;
140  const EntityHandle *conn;
141  CHKERR m_field.get_moab().get_connectivity(*it, conn, num_nodes, true);
142  int order[] = {0, 1, 2};
143  Range adj_tets;
144  CHKERR m_field.get_moab().get_adjacencies(&*it, 1, 3, true, adj_tets);
145  adj_tets = intersect(adj_tets, tets);
146  if (adj_tets.size() == 1) {
147  int side_number;
148  int sense;
149  int offset;
150  CHKERR m_field.get_moab().side_number(adj_tets[0], *it, side_number,
151  sense, offset);
152  if (sense == -1) {
153  order[0] = 1;
154  order[1] = 0;
155  }
156  }
157  tetgen_moab_map[MBTRI | (ii << MB_TYPE_WIDTH)] = *it;
158  moab_tetgen_map[*it] = MBTRI | (ii << MB_TYPE_WIDTH);
159  for (int nn = 0; nn < 3; nn++) {
160  if (moab_tetgen_map.find(conn[order[nn]]) == moab_tetgen_map.end()) {
161  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
162  "data inconsistency between TetGen and MoAB");
163  }
164  in.trifacelist[3 * ii + nn] =
165  moab_tetgen_map[conn[order[nn]]] >> MB_TYPE_WIDTH;
166  }
167  }
168  }
169 
170  Range edges = ents.subset_by_type(MBEDGE);
171  in.numberofedges = edges.size();
172  if (in.numberofedges > 0) {
173  in.edgelist = new int[2 * in.numberofedges];
174  in.edgemarkerlist = new int[in.numberofedges];
175  // std::fill(&in.edgemarkerlist[0],&in.edgemarkerlist[in.numberofedges],1);
176  CHKERR m_field.get_moab().tag_get_data(th_marker, edges, in.edgemarkerlist);
177  it = edges.begin();
178  for (int ii = 0; it != edges.end(); it++, ii++) {
179  int num_nodes;
180  const EntityHandle *conn;
181  CHKERR m_field.get_moab().get_connectivity(*it, conn, num_nodes, true);
182  tetgen_moab_map[MBEDGE | (ii << MB_TYPE_WIDTH)] = *it;
183  moab_tetgen_map[*it] = MBEDGE | (ii << MB_TYPE_WIDTH);
184  for (int nn = 0; nn < 2; nn++) {
185  if (moab_tetgen_map.find(conn[nn]) == moab_tetgen_map.end()) {
186  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
187  "data inconsistency between TetGen and MoAB");
188  }
189  in.edgelist[2 * ii + nn] = moab_tetgen_map[conn[nn]] >> MB_TYPE_WIDTH;
190  }
191  }
192  }
193 
195 }
196 
198  moabTetGen_Map &moab_tetgen_map,
199  tetGenMoab_Map &tetgen_moab_map,
200  std::map<int, Range> &type_ents) {
202 
203  Interface &m_field = cOre;
204  //
205  // ErrorCode rval;
206  in.pointparamlist = new tetgenio::pointparam[in.numberofpoints];
207  // std::vector<bool> points_is_set(in.numberofpoints,false);
208  std::map<int, Range>::iterator mit = type_ents.begin();
209  for (; mit != type_ents.end(); mit++) {
210  if (mit->first < 0 && mit->first > 3) {
211  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
212  "Wrong TetGen point type");
213  }
214  Range::iterator it = mit->second.begin();
215  for (; it != mit->second.end(); it++) {
216  moabTetGen_Map::iterator miit = moab_tetgen_map.find(*it);
217  if (miit == moab_tetgen_map.end()) {
218  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
219  "Data inconsistency between TetGen and MoAB");
220  continue;
221  }
222  int id = miit->second >> MB_TYPE_WIDTH;
223  in.pointparamlist[id].uv[0] = 0;
224  in.pointparamlist[id].uv[1] = 0;
225  in.pointparamlist[id].type = mit->first;
226  in.pointparamlist[id].tag = m_field.get_moab().id_from_handle(*it) + 1;
227  // points_is_set[id] = true;
228  }
229  }
230 
231  // // Check only if tag and type is set to all points
232  // for(
233  // std::vector<bool>::iterator bit = points_is_set.begin();
234  // bit!=points_is_set.end();bit++
235  // ) {
236  // if(!*bit) {
237  // SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"Point type for
238  // TetGen is not set");
239  // }
240  // }
241 
243 }
244 
246 TetGenInterface::outData(tetgenio &in, tetgenio &out,
247  std::map<EntityHandle, unsigned long> &moab_tetgen_map,
248  std::map<unsigned long, EntityHandle> &tetgen_moab_map,
249  Range *ents, bool id_in_tags, bool error_if_created,
250  bool assume_first_nodes_the_same, Tag th) {
252 
253  Interface &m_field = cOre;
254 
255  Tag th_marker;
256  int def_marker = 0;
257  CHKERR m_field.get_moab().tag_get_handle(
258  "TETGEN_MARKER", 1, MB_TYPE_INTEGER, th_marker,
259  MB_TAG_CREAT | MB_TAG_SPARSE, &def_marker);
260 
261  int num_nodes = 0;
262 
263  int ii = 0;
264  for (; ii < out.numberofpoints; ii++) {
265  if (ii < in.numberofpoints) {
266  bool mem_the_same = memcmp(&in.pointlist[3 * ii], &out.pointlist[3 * ii],
267  3 * sizeof(double)) == 0;
268  if (assume_first_nodes_the_same || mem_the_same) {
269  unsigned long iii = MBVERTEX | (ii << MB_TYPE_WIDTH);
270  if (tetgen_moab_map.find(iii) != tetgen_moab_map.end()) {
271  if (!mem_the_same) {
272  if (th)
273  CHKERR m_field.get_moab().tag_set_data(th, &tetgen_moab_map[iii],
274  1, &out.pointlist[3 * ii]);
275  else
276  CHKERR m_field.get_moab().set_coords(&tetgen_moab_map[iii], 1,
277  &out.pointlist[3 * ii]);
278  }
279  CHKERR m_field.get_moab().tag_set_data(
280  th_marker, &tetgen_moab_map[iii], 1, &out.pointmarkerlist[ii]);
281  continue;
282  } else {
283  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
284  "data inconsistency between TetGen and MoAB");
285  }
286  }
287  }
288  if (id_in_tags) {
289  if (out.pointparamlist[ii].tag > 0) {
290  EntityHandle node;
291  CHKERR m_field.get_moab().handle_from_id(
292  MBVERTEX, in.pointparamlist[ii].tag - 1, node);
293  if (moab_tetgen_map.find(node) != moab_tetgen_map.end()) {
294  continue;
295  }
296  }
297  }
298  if (error_if_created) {
299  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
300  "node should not be created");
301  }
302  num_nodes++;
303  }
304 
305  ReadUtilIface *iface;
306  CHKERR m_field.get_moab().query_interface(iface);
307 
308  if (num_nodes) {
309  vector<double *> arrays;
310  EntityHandle startv;
311  CHKERR iface->get_node_coords(3, num_nodes, 0, startv, arrays);
312  Range verts(startv, startv + num_nodes - 1);
313  int ii = in.numberofpoints;
314  for (Range::iterator vit = verts.begin(); vit != verts.end(); vit++, ii++) {
315  for (int nn = 0; nn != 3; nn++)
316  arrays[nn][ii - in.numberofpoints] = out.pointlist[3 * ii + nn];
317  moab_tetgen_map[*vit] = MBVERTEX | (ii << MB_TYPE_WIDTH);
318  tetgen_moab_map[MBVERTEX | (ii << MB_TYPE_WIDTH)] = *vit;
319  if (ents != NULL)
320  ents->insert(*vit);
321  }
322  CHKERR m_field.get_moab().tag_set_data(
323  th_marker, verts, &out.pointmarkerlist[in.numberofpoints]);
324  if(th)
325  CHKERR m_field.get_moab().tag_set_data(
326  th, verts, &out.pointlist[3 * in.numberofpoints]);
327  }
328 
329  std::vector<int> tetgen_ii;
330 
331  // Build tets
332  std::vector<EntityHandle> conn_seq_tet;
333  conn_seq_tet.reserve(4 * out.numberoftetrahedra);
334  tetgen_ii.reserve(out.numberoftetrahedra);
335  conn_seq_tet.clear();
336  tetgen_ii.clear();
337  ii = 0;
338  for (; ii < out.numberoftetrahedra; ii++) {
339  unsigned long iii = MBTET | (ii << MB_TYPE_WIDTH);
340  if (ii < in.numberoftetrahedra) {
341  if (memcmp(&in.tetrahedronlist[4 * ii], &out.tetrahedronlist[4 * ii],
342  4 * sizeof(int)) == 0) {
343  if (tetgen_moab_map.find(iii) != tetgen_moab_map.end()) {
344  if (ents != NULL)
345  ents->insert(tetgen_moab_map[iii]);
346  continue;
347  } else {
348  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
349  "data inconsistency between TetGen and MoAB");
350  }
351  }
352  }
353  EntityHandle conn[4];
354  for (int nn = 0; nn < 4; nn++) {
355  int nnn = out.tetrahedronlist[4 * ii + nn];
356  if (tetgen_moab_map.find(MBVERTEX | (nnn << MB_TYPE_WIDTH)) ==
357  tetgen_moab_map.end()) {
358  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
359  "data inconsistency between TetGen and MoAB");
360  }
361  conn[nn] = tetgen_moab_map.at(MBVERTEX | (nnn << MB_TYPE_WIDTH));
362  }
363  // auto get_coords = [&m_field](const EntityHandle *conn) {
364  // VectorDouble12 coords(12);
365  // CHKERR m_field.get_moab().get_coords(conn, 4, &coords[0]);
366  // return coords;
367  // };
368  // auto get_volume = [](VectorDouble12 &coords) {
369  // return Tools::tetVolume(&coords[0]);
370  // };
371  // if (get_volume(get_coords(conn)) < 0) {
372  // EntityHandle n0 = conn[0];
373  // conn[0] = conn[1];
374  // conn[1] = n0;
375  // }
376  Range tets;
377  CHKERR m_field.get_moab().get_adjacencies(conn, 4, 3, false, tets);
378  bool tet_found = false;
379  for (auto tet : tets) {
380  const EntityHandle *tet_conn;
381  int num_nodes;
382  CHKERR m_field.get_moab().get_connectivity(tet, tet_conn, num_nodes,
383  true);
384  const EntityHandle *p = std::find(tet_conn, &tet_conn[4], conn[0]);
385  if (p != &tet_conn[4]) {
386  int s = std::distance(tet_conn, p);
387  int n = 0;
388  for (; n != 4; ++n) {
389  const int idx[] = {0, 1, 2, 3, 0, 1, 2, 3};
390  if (tet_conn[idx[s + n]] != conn[n])
391  break;
392  }
393  if (n == 4 && !tet_found) {
394  moab_tetgen_map[tet] = iii;
395  tetgen_moab_map[iii] = tet;
396  if (ents != NULL)
397  ents->insert(tet);
398  tet_found = true;
399  } else if (n == 4) {
400  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
401  "More that one tet with the same connectivity");
402  }
403  } else {
404  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "Impossible case");
405  }
406  }
407  if (!tet_found) {
408  for (int nn = 0; nn != 4; nn++) {
409  conn_seq_tet.push_back(conn[nn]);
410  }
411  tetgen_ii.push_back(ii);
412  }
413  }
414 
415  Range new_tets;
416  if (!conn_seq_tet.empty()) {
417  EntityHandle starte; // Connectivity
418  EntityHandle *conn;
419  int num_el = tetgen_ii.size();
420  CHKERR iface->get_element_connect(num_el, 4, MBTET, 0, starte, conn);
421  std::copy(conn_seq_tet.begin(), conn_seq_tet.end(), conn);
422  CHKERR iface->update_adjacencies(starte, num_el, 4, conn);
423  new_tets = Range(starte, starte + num_el - 1);
424  std::vector<int>::iterator ii_it = tetgen_ii.begin();
425  int ii = 0;
426  for (Range::iterator tit = new_tets.begin(); tit != new_tets.end();
427  tit++, ii_it++, ii++) {
428  unsigned long iii = MBTET | ((*ii_it) << MB_TYPE_WIDTH);
429  moab_tetgen_map[*tit] = iii;
430  tetgen_moab_map[iii] = *tit;
431  }
432  if (ents != NULL)
433  ents->merge(new_tets);
434  }
435 
437 }
439 TetGenInterface::outData(tetgenio &in, tetgenio &out,
440  std::map<EntityHandle, unsigned long> &moab_tetgen_map,
441  std::map<unsigned long, EntityHandle> &tetgen_moab_map,
442  BitRefLevel bit, bool id_in_tags,
443  bool error_if_created,
444  bool assume_first_nodes_the_same, Tag th) {
446  Interface &m_field = cOre;
447  Range ents;
448  CHKERR outData(in, out, moab_tetgen_map, tetgen_moab_map, &ents, id_in_tags,
449  error_if_created, assume_first_nodes_the_same,th);
450  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
451  ents.subset_by_type(MBTET), bit);
453 }
455  std::vector<std::pair<Range, int>> &markers, tetgenio &in,
456  std::map<EntityHandle, unsigned long> &moab_tetgen_map,
457  std::map<unsigned long, EntityHandle> &tetgen_moab_map) {
459  ErrorCode rval;
460  Interface &m_field = cOre;
461  in.numberoffacets = markers.size();
462  in.facetlist = new tetgenio::facet[in.numberoffacets];
463  in.facetmarkerlist = new int[in.numberoffacets];
464  std::vector<std::pair<Range, int>>::iterator mit = markers.begin();
465  for (int ii = 0; mit != markers.end(); mit++, ii++) {
466  in.facetmarkerlist[ii] = mit->second;
467  Range &faces = mit->first;
468  tetgenio::facet *f = &(in.facetlist[ii]);
469  f->numberofpolygons = faces.size();
470  f->polygonlist = new tetgenio::polygon[f->numberofpolygons];
471  int jj = 0;
472  for (int dd = 3; dd >= 0; dd--) {
473  Range dd_faces = faces.subset_by_dimension(dd);
474  if (dd_faces.empty())
475  continue;
476  Range::iterator it = dd_faces.begin();
477  for (; it != dd_faces.end(); it++, jj++) {
478  int num_nodes;
479  const EntityHandle *conn;
480  tetgenio::polygon *p = &(f->polygonlist[jj]);
481  switch (type_from_handle(*it)) {
482  case MBVERTEX: {
483  p->numberofvertices = 1;
484  conn = &*it;
485  } break;
486  default: {
487  rval =
488  m_field.get_moab().get_connectivity(*it, conn, num_nodes, true);
490  p->numberofvertices = num_nodes;
491  }
492  }
493  p->vertexlist = new int[p->numberofvertices];
494  for (int nn = 0; nn < p->numberofvertices; nn++) {
495  if (moab_tetgen_map.find(conn[nn]) == moab_tetgen_map.end()) {
496  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
497  "data inconsistency between TetGen and MoAB");
498  }
499  p->vertexlist[nn] = moab_tetgen_map[conn[nn]] >> MB_TYPE_WIDTH;
500  }
501  }
502  }
503  // holes
504  f->numberofholes = 0;
505  f->holelist = NULL;
506  }
508 }
510  std::map<EntityHandle, unsigned long> &tetgen_moab_map, tetgenio &out,
511  Range *ents, idxRange_Map *ents_map, bool only_non_zero) {
513  ErrorCode rval;
514  Interface &m_field = cOre;
515  Tag th_marker;
516  int def_marker = 0;
517  rval = m_field.get_moab().tag_get_handle(
518  "TETGEN_MARKER", 1, MB_TYPE_INTEGER, th_marker,
519  MB_TAG_CREAT | MB_TAG_SPARSE, &def_marker);
521  int ii = 0;
522  for (; ii < out.numberoftrifaces; ii++) {
523  if (only_non_zero) {
524  if (out.trifacemarkerlist[ii] == 0) {
525  continue;
526  }
527  }
528  EntityHandle conn[3];
529  for (int nn = 0; nn < 3; nn++) {
530  int iii = MBVERTEX | (out.trifacelist[3 * ii + nn] << MB_TYPE_WIDTH);
531  if (tetgen_moab_map.find(iii) == tetgen_moab_map.end()) {
532  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
533  "data inconsistency between TetGen and MoAB");
534  } else {
535  conn[nn] = tetgen_moab_map[iii];
536  }
537  }
538  Range face;
539  rval = m_field.get_moab().get_adjacencies(conn, 3, 2, false, face);
541  face = face.subset_by_type(MBTRI);
542  if (face.size() != 1) {
543  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
544  "data inconsistency between TetGen and MoAB, %u", face.size());
545  }
546  if (ents != NULL)
547  ents->merge(face);
548  if (ents_map != NULL)
549  (*ents_map)[out.trifacemarkerlist[ii]].merge(face);
550  rval = m_field.get_moab().tag_set_data(th_marker, &*face.begin(), 1,
551  &out.trifacemarkerlist[ii]);
553  }
555 }
557  std::vector<std::pair<EntityHandle, int>> &regions, tetgenio &in, Tag th) {
559  ErrorCode rval;
560  Interface &m_field = cOre;
561  in.numberofregions = regions.size();
562  in.regionlist = new double[5 * in.numberofregions];
563  int kk = 0;
564  std::vector<std::pair<EntityHandle, int>>::iterator it = regions.begin();
565  for (int ii = 0; it != regions.end(); it++, ii++) {
566  double coords[3];
567  switch (type_from_handle(it->first)) {
568  case MBVERTEX: {
569  if (th) {
570  rval = m_field.get_moab().tag_get_data(th, &it->first, 1, coords);
572  } else {
573  rval = m_field.get_moab().get_coords(&it->first, 1, coords);
575  }
576  } break;
577  case MBTET: {
578  int num_nodes;
579  const EntityHandle *conn;
580  rval =
581  m_field.get_moab().get_connectivity(it->first, conn, num_nodes, true);
583  double _coords[12];
584  if (th) {
585  rval = m_field.get_moab().tag_get_data(th, conn, num_nodes, _coords);
587  } else {
588  rval = m_field.get_moab().get_coords(conn, num_nodes, _coords);
590  }
591  coords[0] = (_coords[0] + _coords[3] + _coords[6] + _coords[9]) / 4.;
592  coords[1] = (_coords[1] + _coords[4] + _coords[7] + _coords[10]) / 4.;
593  coords[2] = (_coords[2] + _coords[5] + _coords[8] + _coords[11]) / 4.;
594  } break;
595  default:
596  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
597  }
598  for (int nn = 0; nn < 3; nn++) {
599  in.regionlist[kk++] = coords[nn];
600  }
601  in.regionlist[kk++] = it->second;
602  in.regionlist[kk++] = it->second;
603  }
605 }
607  std::map<EntityHandle, unsigned long> &tetgen_moab_map, tetgenio &out,
608  Range *ents, idxRange_Map *ents_map) {
610  Interface &m_field = cOre;
611  int nbattributes = out.numberoftetrahedronattributes;
612  if (nbattributes == 0) {
613  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
614  "tetgen has no regions attributes");
615  }
616  Tag th_region;
617  rval = m_field.get_moab().tag_get_handle("TETGEN_REGION", th_region);
618  if (rval == MB_SUCCESS) {
619  rval = m_field.get_moab().tag_delete(th_region);
621  }
622  double def_marker = 0;
623  CHKERR m_field.get_moab().tag_get_handle(
624  "TETGEN_REGION", nbattributes, MB_TYPE_DOUBLE, th_region,
625  MB_TAG_CREAT | MB_TAG_SPARSE, &def_marker);
626  int ii = 0;
627  for (; ii < out.numberoftetrahedra; ii++) {
628  int jj = 0;
629  for (; jj < nbattributes; jj++) {
630  double id = out.tetrahedronattributelist[ii * nbattributes + jj];
631  int iii = MBTET | (ii << MB_TYPE_WIDTH);
632  if (tetgen_moab_map.find(iii) == tetgen_moab_map.end()) {
633  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
634  "data inconsistency between TetGen and MoAB");
635  }
636  EntityHandle ent = tetgen_moab_map[iii];
637  CHKERR m_field.get_moab().tag_set_data(th_region, &ent, 1, &id);
638  if (ents != NULL)
639  ents->insert(ent);
640  if (ents_map != NULL)
641  (*ents_map)[id].insert(ent);
642  }
643  }
645 }
646 MoFEMErrorCode TetGenInterface::tetRahedralize(char switches[], tetgenio &in,
647  tetgenio &out) {
649  tetgenbehavior a;
650  a.parse_commandline(switches);
651  tetrahedralize(&a, &in, &out);
653 }
654 MoFEMErrorCode TetGenInterface::loadPoly(char file_name[], tetgenio &in) {
656  if (!in.load_poly(file_name)) {
657  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
658  "can not read TetGen poly file");
659  }
661 }
663  bool *result,
664  const double eps) {
666  double *pa = &coords[0];
667  double *pb = &coords[3];
668  double *pc = &coords[6];
669  double *pd = &coords[9];
670  double adx = pa[0] - pd[0];
671  double bdx = pb[0] - pd[0];
672  double cdx = pc[0] - pd[0];
673  double ady = pa[1] - pd[1];
674  double bdy = pb[1] - pd[1];
675  double cdy = pc[1] - pd[1];
676  double adz = pa[2] - pd[2];
677  double bdz = pb[2] - pd[2];
678  double cdz = pc[2] - pd[2];
679  double v = adx * (bdy * cdz - bdz * cdy) + bdx * (cdy * adz - cdz * ady) +
680  cdx * (ady * bdz - adz * bdy);
681  double l = sqrt(pow(pa[0] - pb[0], 2) + pow(pa[1] - pb[1], 2) +
682  pow(pa[2] - pb[2], 2)) +
683  sqrt(pow(pa[0] - pc[0], 2) + pow(pa[1] - pc[1], 2) +
684  pow(pa[2] - pc[2], 2)) +
685  sqrt(pow(pa[0] - pd[0], 2) + pow(pa[1] - pd[1], 2) +
686  pow(pa[2] - pd[2], 2)) +
687  sqrt(pow(pb[0] - pc[0], 2) + pow(pb[1] - pc[1], 2) +
688  pow(pb[2] - pc[2], 2)) +
689  sqrt(pow(pb[0] - pd[0], 2) + pow(pb[1] - pd[1], 2) +
690  pow(pb[2] - pd[2], 2)) +
691  sqrt(pow(pc[0] - pd[0], 2) + pow(pc[1] - pd[1], 2) +
692  pow(pc[2] - pd[2], 2));
693  // std::cerr << fabs(v/pow(l,3)) << " ";
694  *result = fabs(v / pow(l, 3)) < eps ? true : false;
696 }
698  std::vector<Range> &sorted,
699  const double eps, Tag th) {
701 
702  Interface &m_field = cOre;
703  Skinner skin(&m_field.get_moab());
704 
705  for (;;) {
706 
707  Range noplanar_to_anyother;
708  std::vector<Range>::iterator vit = sorted.begin();
709 
710  do {
711 
712  bool repeat = false;
713 
714  // get edges on vit skin
715  Range skin_edges;
716  CHKERR skin.find_skin(0, *vit, false, skin_edges);
717 
718  // skin edge nodes
719  Range skin_edges_nodes;
720  CHKERR m_field.get_moab().get_connectivity(skin_edges, skin_edges_nodes,
721  true);
722 
723  // get tris adjacent to vit skin edges
724  Range skin_edges_tris;
725  CHKERR m_field.get_moab().get_adjacencies(
726  skin_edges, 2, false, skin_edges_tris, moab::Interface::UNION);
727  // get tris which are part of facet
728  Range inner_tris = intersect(skin_edges_tris, *vit);
729  Range outer_tris = intersect(skin_edges_tris, tris);
730 
731  // tris coplanar with vit tris
732  Range coplanar_tris;
733 
734  Range::iterator tit = outer_tris.begin();
735  for (; tit != outer_tris.end(); tit++) {
736  Range tit_conn;
737  CHKERR m_field.get_moab().get_connectivity(&*tit, 1, tit_conn, true);
738  tit_conn = subtract(tit_conn, skin_edges_nodes);
739  if (tit_conn.empty()) {
740  coplanar_tris.insert(*tit);
741  noplanar_to_anyother.erase(*tit);
742  repeat = true;
743  } else {
744  Range tit_edges;
745  CHKERR m_field.get_moab().get_adjacencies(
746  &*tit, 1, 1, false, tit_edges, moab::Interface::UNION);
747  tit_edges = intersect(tit_edges, skin_edges);
748  if (tit_edges.size() != 1) {
749  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
750  "data inconsistency");
751  }
752  Range inner_tri;
753  CHKERR m_field.get_moab().get_adjacencies(
754  tit_edges, 2, false, inner_tri, moab::Interface::UNION);
755  inner_tri = intersect(inner_tri, inner_tris);
756  if (inner_tri.size() != 1) {
757  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
758  "data inconsistency");
759  }
760  // get connectivity
761  int num_nodes;
762  const EntityHandle *inner_tri_conn;
763  CHKERR m_field.get_moab().get_connectivity(
764  *inner_tri.begin(), inner_tri_conn, num_nodes, true);
765  double coords[12];
766  if (th) {
767  CHKERR m_field.get_moab().tag_get_data(th, inner_tri_conn, 3,
768  coords);
769  CHKERR m_field.get_moab().tag_get_data(th, &*tit_conn.begin(), 1,
770  &coords[9]);
771  } else {
772  CHKERR m_field.get_moab().get_coords(inner_tri_conn, 3, coords);
773  CHKERR m_field.get_moab().get_coords(&*tit_conn.begin(), 1,
774  &coords[9]);
775  }
776  bool coplanar;
777  CHKERR checkPlanar_Trinagle(coords, &coplanar, eps);
778  if (coplanar) {
779  coplanar_tris.insert(*tit);
780  noplanar_to_anyother.erase(*tit);
781  repeat = true;
782  } else {
783  noplanar_to_anyother.insert(*tit);
784  }
785  }
786  }
787 
788  vit->merge(coplanar_tris);
789  tris = subtract(tris, *vit);
790 
791  if (repeat) {
792  vit = sorted.begin();
793  } else {
794  vit++;
795  }
796 
797  } while (vit != sorted.end());
798 
799  if (noplanar_to_anyother.empty()) {
801  } else {
802  Range seed;
803  seed.insert(noplanar_to_anyother[0]);
804  tris.erase(noplanar_to_anyother[0]);
805  sorted.push_back(seed);
806  }
807  }
808 
810 }
811 
813  Range &tris, std::vector<std::vector<Range>> &sorted, const double eps) {
815 
816  // PetscAttachDebugger();
817 
818  Range seed;
819  seed.insert(tris[0]);
820  tris.erase(tris[0]);
821  std::vector<Range> vec_seed;
822  vec_seed.push_back(seed);
823  sorted.push_back(vec_seed);
824 
825  for (;;) {
826  std::vector<Range> &vec = sorted.back();
827  CHKERR groupPlanar_Triangle(tris, vec, eps);
828  if (tris.empty()) {
830  } else {
831  Range seed;
832  seed.insert(tris[0]);
833  tris.erase(tris[0]);
834  std::vector<Range> vec_seed;
835  vec_seed.push_back(seed);
836  sorted.push_back(vec_seed);
837  }
838  }
839 
841 }
843  bool reduce_edges,
844  Range *not_reducable_nodes,
845  const double eps, Tag th) {
847  // FIXME: assumes that are no holes
848 
849  if (ents.empty()) {
850  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
851  "no ents to build polygon");
852  }
853 
854  Interface &m_field = cOre;
855 
856  Skinner skin(&m_field.get_moab());
857 
858  Range skin_edges;
859  CHKERR skin.find_skin(0, ents, false, skin_edges);
860 
861  std::vector<EntityHandle> polygon_nodes;
862  EntityHandle seed = skin_edges[0];
863  Range seen_edges;
864  seen_edges.insert(seed);
865  skin_edges.erase(seed);
866  int num_nodes;
867  const EntityHandle *conn;
868  CHKERR m_field.get_moab().get_connectivity(seed, conn, num_nodes, true);
869  polygon_nodes.push_back(conn[0]);
870  // std::cerr << std::endl;
871  // std::cerr << conn[0] << " " << conn[1] << std::endl;
872  do {
873  EntityHandle last_node = polygon_nodes.back();
874  Range adj_edges;
875  CHKERR m_field.get_moab().get_adjacencies(&last_node, 1, 1, false,
876  adj_edges);
877  adj_edges = intersect(adj_edges, skin_edges);
878  if (adj_edges.size() == 0) {
879  break;
880  }
881  if (adj_edges.size() != 1) {
882  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
883  "should be only one edge");
884  }
885  seen_edges.insert(adj_edges[0]);
886  skin_edges.erase(adj_edges[0]);
887  CHKERR m_field.get_moab().get_connectivity(adj_edges[0], conn, num_nodes,
888  true);
889  EntityHandle add_node = (last_node == conn[0]) ? conn[1] : conn[0];
890  polygon_nodes.push_back(add_node);
891  // std::cerr << "\t" << add_node << std::endl;
892  } while (1);
893 
894  if (reduce_edges) {
895  // std::cerr << "polygon " << polygon_nodes.size();
896  std::vector<EntityHandle>::iterator pit = polygon_nodes.begin();
897  // std::cerr << std::endl;
898  for (; pit != polygon_nodes.end();) {
899  if (not_reducable_nodes != NULL) {
900  if (not_reducable_nodes->find(*pit) != not_reducable_nodes->end()) {
901  pit++;
902  continue;
903  }
904  }
905  EntityHandle mm;
906  if (pit == polygon_nodes.begin()) {
907  mm = polygon_nodes.back();
908  } else {
909  mm = *(pit - 1);
910  }
911  EntityHandle mc = *pit;
912  EntityHandle mp;
913  if (polygon_nodes.back() == mc) {
914  mp = polygon_nodes[0];
915  } else {
916  mp = *(pit + 1);
917  }
918  double coords[9];
919  if (th) {
920  CHKERR m_field.get_moab().tag_get_data(th, &mm, 1, &coords[3 * 0]);
921  CHKERR m_field.get_moab().tag_get_data(th, &mc, 1, &coords[3 * 1]);
922  CHKERR m_field.get_moab().tag_get_data(th, &mp, 1, &coords[3 * 2]);
923  } else {
924  CHKERR m_field.get_moab().get_coords(&mm, 1, &coords[3 * 0]);
925  CHKERR m_field.get_moab().get_coords(&mc, 1, &coords[3 * 1]);
926  CHKERR m_field.get_moab().get_coords(&mp, 1, &coords[3 * 2]);
927  }
928  cblas_daxpy(3, -1, &coords[3 * 1], 1, &coords[3 * 0], 1); // mm = mm -
929  // mc
930  cblas_daxpy(3, -1, &coords[3 * 1], 1, &coords[3 * 2], 1); // mp = mp -
931  // mc
932  double spin[9];
933  CHKERR Spin(spin, &coords[3 * 0]);
934  double l0 = cblas_dnrm2(3, &coords[3 * 0], 1);
935  cblas_dgemv(CblasRowMajor, CblasNoTrans, 3, 3, 1. / l0, spin, 3,
936  &coords[3 * 2], 1, 0., &coords[3 * 1], 1);
937  double dot = cblas_dnrm2(3, &coords[3 * 1], 1);
938  // std::cerr << mm << " " << mc << " " << mp << " " << dot << std::endl;
939  if (dot < eps) {
940  polygon_nodes.erase(pit);
941  pit = polygon_nodes.begin();
942  // std::cerr << std::endl;
943  } else {
944  pit++;
945  }
946  }
947  }
948 
949  Range existing_polygon;
950  CHKERR m_field.get_moab().get_adjacencies(
951  &polygon_nodes[0], polygon_nodes.size(), 2, true, existing_polygon);
952  if (existing_polygon.empty()) {
953  EntityHandle polygon;
954  CHKERR m_field.get_moab().create_element(MBPOLYGON, &polygon_nodes[0],
955  polygon_nodes.size(), polygon);
956  polygons.insert(polygon);
957  } else {
958  polygons.merge(existing_polygon);
959  }
960 
962 }
963 } // namespace MoFEM
964 
965 #endif // TETGEN
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::TetGenInterface::cOre
MoFEM::Core & cOre
Definition: TetGenInterface.hpp:28
MoFEM::TetGenInterface::tetGenMoab_Map
std::map< unsigned long, EntityHandle > tetGenMoab_Map
Definition: TetGenInterface.hpp:33
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
MoFEM::TetGenInterface::getRegionData
MoFEMErrorCode getRegionData(tetGenMoab_Map &tetgen_moab_map, tetgenio &out, Range *ents=NULL, idxRange_Map *ents_map=NULL)
get region data to tetrahedral
Definition: TetGenInterface.cpp:606
CHKERRQ_MOAB
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:454
MB_TYPE_WIDTH
#define MB_TYPE_WIDTH
Definition: definitions.h:226
MoFEM::TetGenInterface::setFaceData
MoFEMErrorCode setFaceData(std::vector< std::pair< Range, int > > &markers, tetgenio &in, moabTetGen_Map &moab_tetgen_map, tetGenMoab_Map &tetgen_moab_map)
set markers to faces
Definition: TetGenInterface.cpp:454
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::TetGenInterface::query_interface
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Definition: TetGenInterface.cpp:31
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
MoFEM::TetGenInterface::moabTetGen_Map
std::map< EntityHandle, unsigned long > moabTetGen_Map
Definition: TetGenInterface.hpp:32
MoFEM::TetGenInterface::tetRahedralize
MoFEMErrorCode tetRahedralize(char switches[], tetgenio &in, tetgenio &out)
run tetgen
Definition: TetGenInterface.cpp:646
MoFEM::TetGenInterface::idxRange_Map
std::map< int, Range > idxRange_Map
Definition: TetGenInterface.hpp:34
FTensor::Tensor2< double, 3, 3 >
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::TetGenInterface::outData
MoFEMErrorCode outData(tetgenio &in, tetgenio &out, moabTetGen_Map &moab_tetgen_map, tetGenMoab_Map &tetgen_moab_map, Range *ents=NULL, bool id_in_tags=false, bool error_if_created=false, bool assume_first_nodes_the_same=false, Tag th=nullptr)
get entities for TetGen data structure
Definition: TetGenInterface.cpp:246
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
a
constexpr double a
Definition: approx_sphere.cpp:30
Spin
PetscErrorCode Spin(double *spinOmega, double *vecOmega)
calculate spin matrix from vector
Definition: fem_tools.c:546
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
MoFEM::TetGenInterface::checkPlanar_Trinagle
MoFEMErrorCode checkPlanar_Trinagle(double coords[], bool *result, const double eps=1e-9)
Definition: TetGenInterface.cpp:662
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
MoFEM::type_from_handle
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1869
t
constexpr double t
plate stiffness
Definition: plate.cpp:59
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
FTensor::Index< 'i', 3 >
convert.n
n
Definition: convert.py:82
MoFEM::TetGenInterface
TetGen interface.
Definition: TetGenInterface.hpp:23
MoFEM::UnknownInterface
base class for all interface classes
Definition: UnknownInterface.hpp:34
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1511
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
Range
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
MoFEM::TetGenInterface::groupRegion_Triangle
MoFEMErrorCode groupRegion_Triangle(Range &tris, std::vector< std::vector< Range > > &sorted, const double eps=1e-9)
Group surface triangles in planar regions.
Definition: TetGenInterface.cpp:812
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
MoFEM::TetGenInterface::groupPlanar_Triangle
MoFEMErrorCode groupPlanar_Triangle(Range &tris, std::vector< Range > &sorted, const double eps=1e-9, Tag th=NULL)
Definition: TetGenInterface.cpp:697
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
ShapeDiffMBTET
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:319
MoFEM::TetGenInterface::inData
MoFEMErrorCode inData(Range &ents, tetgenio &in, moabTetGen_Map &moab_tetgen_map, tetGenMoab_Map &tetgen_moab_map, Tag th=NULL)
create TetGen data structure form range of moab entities
Definition: TetGenInterface.cpp:38
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
FTensor.hpp
Tensors class implemented by Walter Landry.
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
MoFEM::TetGenInterface::setRegionData
MoFEMErrorCode setRegionData(std::vector< std::pair< EntityHandle, int > > &regions, tetgenio &in, Tag th=NULL)
set region data to tetrahedral
Definition: TetGenInterface.cpp:556
MoFEM::TetGenInterface::loadPoly
MoFEMErrorCode loadPoly(char file_name[], tetgenio &in)
load poly file
Definition: TetGenInterface.cpp:654
MoFEM::TetGenInterface::setGeomData
MoFEMErrorCode setGeomData(tetgenio &in, moabTetGen_Map &moab_tetgen_map, tetGenMoab_Map &tetgen_moab_map, std::map< int, Range > &type_ents)
set point tags and type
Definition: TetGenInterface.cpp:197
MoFEM::TetGenInterface::getTriangleMarkers
MoFEMErrorCode getTriangleMarkers(tetGenMoab_Map &tetgen_moab_map, tetgenio &out, Range *ents=NULL, idxRange_Map *ents_map=NULL, bool only_non_zero=true)
get markers to faces
Definition: TetGenInterface.cpp:509
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MoFEM::TetGenInterface::makePolygonFacet
MoFEMErrorCode makePolygonFacet(Range &ents, Range &polygons, bool reduce_edges=false, Range *not_reducable_nodes=NULL, const double eps=1e-9, Tag th=NULL)
Definition: TetGenInterface.cpp:842
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21