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