v0.14.0
NodeMerger.cpp
Go to the documentation of this file.
1 /** \file NodeMerger.cpp
2  * \brief Interface for merging nodes
3  */
4 
5 
6 
7 namespace MoFEM {
8 
10 NodeMergerInterface::query_interface(boost::typeindex::type_index type_index,
11  UnknownInterface **iface) const {
12  *iface = const_cast<NodeMergerInterface *>(this);
13  return 0;
14 }
15 
16 static auto min_non_abs(const double a, const double b) {
17  return std::min(a, b);
18 };
19 
21  : cOre(const_cast<MoFEM::Core &>(core)), minQualityFunction(min_non_abs),
22  successMerge(false), errorIfNoCommonEdge(false) {}
23 
25  Interface &m_field = cOre;
27  ierr = PetscOptionsBegin(m_field.get_comm(), "node_merge",
28  "Node merge options", "none");
29  CHKERRQ(ierr);
30 
31  ierr = PetscOptionsEnd();
32  CHKERRQ(ierr);
33 
35 }
36 
38  EntityHandle father, EntityHandle mother, Range &out_tets, Range *tets_ptr,
39  const bool only_if_improve_quality, const double move,
40  const int line_search, Tag th, int verb) {
41  Interface &m_field = cOre;
44 
45  // Get edges adjacent to father and mother, i.e. mother is merged to father.
46  Range father_edges;
47  CHKERR m_field.get_moab().get_adjacencies(&father, 1, 1, false, father_edges);
48  Range mother_edges;
49  CHKERR m_field.get_moab().get_adjacencies(&mother, 1, 1, false, mother_edges);
50 
51  // Get tets adjacent to mother and father
52  Range father_tets;
53  CHKERR m_field.get_moab().get_adjacencies(&father, 1, 3, false, father_tets);
54  Range mother_tets;
55  CHKERR m_field.get_moab().get_adjacencies(&mother, 1, 3, false, mother_tets);
56  if (tets_ptr != NULL) {
57  mother_tets = intersect(mother_tets, *tets_ptr);
58  father_tets = intersect(father_tets, *tets_ptr);
59  }
60 
61  // Find common edge
62  Range common_edge;
63  common_edge = intersect(father_edges, mother_edges);
64  if (tets_ptr != NULL) {
65  Range tets = unite(father_tets, mother_tets);
66  Range tets_edges;
67  CHKERR m_field.get_moab().get_adjacencies(tets, 1, false, tets_edges,
68  moab::Interface::UNION);
69  common_edge = intersect(common_edge, tets_edges);
70  father_edges = intersect(father_edges, tets_edges);
71  mother_edges = intersect(mother_edges, tets_edges);
72  }
73 
74  // No common edge, merge no possible
75  if (errorIfNoCommonEdge && common_edge.empty()) {
76  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
77  "no common edge between nodes");
78  } else if (common_edge.empty()) {
79  Range seed_tets;
80  if (tets_ptr != NULL) {
81  seed_tets.merge(*tets_ptr);
82  }
83  out_tets = seed_tets;
84  successMerge = false;
86  }
87 
88  // Common edge tets, that tests will be squashed
89  Range edge_tets;
90  CHKERR m_field.get_moab().get_adjacencies(common_edge, 3, true, edge_tets);
91  // Intersect with ptr_tets (usually associated with some bit level)
92  if (tets_ptr != NULL) {
93  edge_tets = intersect(edge_tets, *tets_ptr);
94  }
95  // Mother tets, has only one mother vertex and no father vertex.
96  mother_tets = subtract(mother_tets, edge_tets);
97  father_tets = subtract(father_tets, edge_tets);
98 
99  auto get_coords = [&m_field](Tag th, const EntityHandle *conn,
100  const int num_nodes) {
101  VectorDouble coords(3 * num_nodes);
102  if (th == NULL) {
103  CHKERR m_field.get_moab().get_coords(conn, num_nodes, &coords[0]);
104  } else {
105  CHKERR m_field.get_moab().tag_get_data(th, conn, num_nodes, &coords[0]);
106  }
107  return coords;
108  };
109 
110  auto get_tensor = [](VectorDouble &coords, const int shift) {
112  &coords[shift], &coords[shift + 1], &coords[shift + 2]);
113  };
114 
115  // move father coord is move > 0
117  if (move > 0) {
118  EntityHandle conn[] = {father, mother};
119  VectorDouble coords = get_coords(th, conn, 2);
120  auto t_n0 = get_tensor(coords, 0);
121  auto t_n1 = get_tensor(coords, 3);
122  t_move(i) = t_n0(i) + move * (t_n1(i) - t_n0(i));
123  }
124 
125  if (line_search > 0) {
126  Range check_tests = unite(father_tets, mother_tets);
127  CHKERR lineSearch(check_tests, father, mother, line_search, t_move, th);
128  }
129 
130  if (only_if_improve_quality) {
131  Range check_tests = mother_tets;
132  // no need to check father tets since no change of quality for them
133  if (move > 0 || line_search) {
134  check_tests.merge(father_tets);
135  }
136 
137  double min_quality0 = 1;
138  CHKERR minQuality(edge_tets, 0, 0, NULL, min_quality0, th,
140  CHKERR minQuality(check_tests, 0, 0, NULL, min_quality0, th,
142  double min_quality = 1;
143  CHKERR minQuality(check_tests, father, mother,
144  ((move > 0) || line_search) ? &t_move(0) : NULL,
145  min_quality, th, minQualityFunction);
146  if (min_quality < min_quality0) {
147  if (tets_ptr != NULL) {
148  out_tets = *tets_ptr;
149  }
150  successMerge = false;
152  }
153  }
154 
155  // Move node
156  if (move > 0 || line_search) {
157  if (th == NULL) {
158  CHKERR m_field.get_moab().set_coords(&father, 1, &t_move(0));
159  } else {
160  CHKERR m_field.get_moab().tag_set_data(th, &father, 1, &t_move(0));
161  }
162  }
163 
164  auto get_conn = [&m_field](const EntityHandle ent,
165  int *ret_num_nodes = nullptr) {
166  int num_nodes;
167  const EntityHandle *conn;
168  CHKERR m_field.get_moab().get_connectivity(ent, conn, num_nodes, true);
169  if (ret_num_nodes)
170  *ret_num_nodes = num_nodes;
171  return conn;
172  };
173 
174  auto create_tet = [this, &m_field](const EntityHandle *new_conn,
175  const EntityHandle parent) {
176  EntityHandle tet;
177  Range tets;
178  CHKERR m_field.get_moab().get_adjacencies(new_conn, 4, 3, false, tets);
179  bool tet_found = false;
180  for (auto it_tet : tets) {
181  const EntityHandle *tet_conn;
182  int num_nodes;
183  CHKERR m_field.get_moab().get_connectivity(it_tet, tet_conn, num_nodes,
184  true);
185  const EntityHandle *p = std::find(tet_conn, &tet_conn[4], new_conn[0]);
186  if (p != tet_conn + 4) {
187  int s = std::distance(tet_conn, p);
188  int n = 0;
189  for (; n != 4; ++n) {
190  const int idx[] = {0, 1, 2, 3, 0, 1, 2, 3};
191  if (tet_conn[idx[s + n]] != new_conn[n])
192  break;
193  }
194  if (n == 4 && !tet_found) {
195  tet = it_tet;
196  tet_found = true;
197  } else if (n == 4) {
198  THROW_MESSAGE("More that one tet with the same connectivity");
199  }
200  }
201  }
202  if (!tet_found) {
203  // Create tet with new connectivity
204  CHKERR m_field.get_moab().create_element(MBTET, new_conn, 4, tet);
205  CHKERR m_field.get_moab().tag_set_data(cOre.get_th_RefParentHandle(),
206  &tet, 1, &parent);
207  parentChildMap.insert(ParentChild(parent, tet));
208  }
209  return tet;
210  };
211 
212  // clear map
213  parentChildMap.clear();
214 
215  Range created_tets;
216  for (auto m_tet : mother_tets) {
217  const EntityHandle *conn = get_conn(m_tet);
218  EntityHandle new_conn[4];
219  // Replace mother vertices by father vertices
220  int nb_mother_verts = 0;
221  for (int nn = 0; nn != 4; ++nn) {
222  if (conn[nn] == father) {
223  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
224  "Tet has father vertex, impossible but here it is");
225  }
226  if (conn[nn] == mother) {
227  new_conn[nn] = father;
228  nb_mother_verts++;
229  } else {
230  new_conn[nn] = conn[nn];
231  }
232  }
233  if (nb_mother_verts != 1) {
234  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
235  "Tet should have only one vertex but have %d", nb_mother_verts);
236  }
237 
238  VectorDouble new_coords = get_coords(th, new_conn, 4);
239 
240  // add tet to range
241  created_tets.insert(create_tet(new_conn, m_tet));
242  }
243 
244  // Loop over mother adjacent entities to use them as parents
245  auto get_adj_ents = [&](const Range &ents, const bool create) {
246  Range adj;
247  for (int dd = 1; dd <= 2; dd++)
248  CHKERR m_field.get_moab().get_adjacencies(ents, dd, create, adj,
249  moab::Interface::UNION);
250  return adj;
251  };
252  auto adj_crated_ents = get_adj_ents(created_tets, true);
253  adj_crated_ents.erase(common_edge[0]);
254 
255  FaceMapIdx face_map;
256  for (auto ent : adj_crated_ents) {
257  int num_nodes;
258  const EntityHandle *conn = get_conn(ent, &num_nodes);
259  EntityHandle small_conn[num_nodes];
260  int ii = 0;
261  bool father_node = false;
262  for (int nn = 0; nn != num_nodes; nn++) {
263  if (conn[nn] == father)
264  father_node = true;
265  else
266  small_conn[ii++] = conn[nn];
267  }
268  if (father_node) {
269  if (ii > 1)
270  std::sort(&small_conn[0], &small_conn[ii]);
271  if (ii == 2)
272  face_map.insert(FaceMap(ent, small_conn[0], small_conn[1]));
273  else
274  face_map.insert(FaceMap(ent, small_conn[0], 0));
275  }
276  }
277 
278  auto adj_mother_ents = get_adj_ents(mother_tets, false);
279  adj_mother_ents.erase(common_edge[0]);
280  for (auto ent : adj_mother_ents) {
281  int num_nodes;
282  const EntityHandle *conn = get_conn(ent, &num_nodes);
283  EntityHandle small_conn[num_nodes];
284  int nb_new_node = 0;
285  int nn = 0;
286  int ii = 0;
287  for (; nn != num_nodes; ++nn) {
288  if (conn[nn] == mother) {
289  nb_new_node++;
290  } else {
291  small_conn[ii++] = conn[nn];
292  }
293  }
294  if (nb_new_node > 0) {
295  if (ii > 1)
296  std::sort(&small_conn[0], &small_conn[ii]);
297 
298  EntityHandle n0 = small_conn[0], n1 = 0;
299  if (ii == 2)
300  n1 = small_conn[1];
301 
302  FaceMapIdx::iterator fit = face_map.find(boost::make_tuple(n0, n1));
303  if (fit != face_map.end()) {
304  const EntityHandle child = fit->e;
305  const EntityHandle parent = ent;
306  if (m_field.get_moab().dimension_from_handle(parent) !=
307  m_field.get_moab().dimension_from_handle(child))
308  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
309  "Huston we have a problem!");
310 
311  // Set parent child relation
312  CHKERR m_field.get_moab().tag_set_data(cOre.get_th_RefParentHandle(),
313  &child, 1, &parent);
314  // Create map
315  parentChildMap.insert(ParentChild(parent, child));
316  }
317  }
318  }
319 
320  // Seed tets to given bit level
321  Range seed_tets;
322  if (tets_ptr != NULL) {
323  seed_tets.merge(*tets_ptr);
324  seed_tets = subtract(seed_tets, edge_tets);
325  seed_tets = subtract(seed_tets, mother_tets);
326  }
327  seed_tets.merge(created_tets);
328  out_tets.swap(seed_tets);
329 
330  successMerge = true;
331 
332  if (verb > VERY_VERBOSE)
333  std::cerr << "Nodes merged" << endl;
334 
336 }
337 
340  EntityHandle mother, double *coords_move,
341  double &min_quality, Tag th,
342  boost::function<double(double, double)> f) {
343  Interface &m_field = cOre;
344  double coords[12];
347  for (Range::iterator tit = check_tests.begin(); tit != check_tests.end();
348  tit++) {
349  const EntityHandle *conn;
350  int num_nodes;
351  CHKERR m_field.get_moab().get_connectivity(*tit, conn, num_nodes, true);
352  if (mother > 0) {
353  EntityHandle new_conn[4];
354  // Replace mother vertices by father vertices
355  int nb_mother_verts = 0;
356  int father_nn = 0;
357  for (int nn = 0; nn < 4; nn++) {
358  if (conn[nn] == father) {
359  father_nn = nn;
360  }
361  if (conn[nn] == mother) {
362  new_conn[nn] = father;
363  father_nn = nn;
364  nb_mother_verts++;
365  } else {
366  new_conn[nn] = conn[nn];
367  }
368  }
369  if (nb_mother_verts > 1) {
370  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
371  "Tet should have no more than one mother vertex but have %d",
372  nb_mother_verts);
373  }
374  if (th == NULL) {
375  CHKERR m_field.get_moab().get_coords(new_conn, num_nodes, coords);
376  } else {
377  CHKERR m_field.get_moab().tag_get_data(th, new_conn, num_nodes, coords);
378  }
379  if (coords_move) {
380  int shift = 3 * father_nn;
381  for (int nn = 0; nn != 3; nn++) {
382  coords[shift + nn] = coords_move[nn];
383  }
384  }
385  } else {
386  if (th == NULL) {
387  CHKERR m_field.get_moab().get_coords(conn, num_nodes, coords);
388  } else {
389  CHKERR m_field.get_moab().tag_get_data(th, conn, num_nodes, coords);
390  }
391  }
392  double quality = Tools::volumeLengthQuality(coords);
393  if (!std::isnormal(quality)) {
394  min_quality = 0;
396  }
397  min_quality = f(min_quality, quality);
398  }
400 };
401 
404  EntityHandle mother, int line_search,
405  FTensor::Tensor1<double, 3> &t_move, Tag th) {
406  Interface &m_field = cOre;
408 
409  EntityHandle conn[] = {father, mother};
410 
411  double coords[6];
412  if (th == NULL) {
413  CHKERR m_field.get_moab().get_coords(conn, 2, coords);
414  } else {
415  CHKERR m_field.get_moab().tag_get_data(th, conn, 2, coords);
416  }
417 
419  FTensor::Tensor1<double, 3> t_coords(coords[0], coords[1], coords[2]);
421  for (int nn = 0; nn != 3; nn++) {
422  t_delta(nn) = coords[3 + nn] - t_coords(nn);
423  }
424 
425  t_move(i) = t_coords(i);
426  double min_quality_i = 1;
427  CHKERR minQuality(check_tests, father, mother, &t_move(0), min_quality_i, th,
429 
430  t_move(i) = t_coords(i) + t_delta(i);
431  double min_quality_k = 1;
432  CHKERR minQuality(check_tests, father, mother, &t_move(0), min_quality_k, th,
434 
435  double alpha_i = 0;
436  double alpha_k = 1;
437 
438  for (int ii = 0; ii != line_search; ii++) {
439  double min_quality = 1;
440  double alpha = (alpha_i + alpha_k) * 0.5;
441  t_move(i) = t_coords(i) + alpha * t_delta(i);
442  CHKERR minQuality(check_tests, father, mother, &t_move(0), min_quality, th,
444  if (min_quality_i >= min_quality_k) {
445  min_quality_k = min_quality;
446  alpha_k = alpha;
447  } else {
448  min_quality_i = min_quality;
449  alpha_i = alpha;
450  }
451  }
452 
453  if (min_quality_i > min_quality_k) {
454  t_move(i) = t_coords(i) + alpha_i * t_delta(i);
455  } else {
456  t_move(i) = t_coords(i) + alpha_k * t_delta(i);
457  }
458 
460 }
461 
463  EntityHandle father, EntityHandle mother, BitRefLevel bit, Range *tets_ptr,
464  const bool only_if_improve_quality, const double move, Tag th) {
465  Interface &m_field = cOre;
467  Range out_tets;
468  CHKERR mergeNodes(father, mother, out_tets, tets_ptr, only_if_improve_quality,
469  move, 0, th);
470  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(out_tets, bit);
472 }
473 
475  EntityHandle father, EntityHandle mother, BitRefLevel bit,
476  BitRefLevel tets_from_bit_ref_level, const bool only_if_improve_quality,
477  const double move, Tag th) {
478  Interface &m_field = cOre;
480  Range level_tets;
481  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
482  tets_from_bit_ref_level, BitRefLevel().set(), MBTET, level_tets);
483  CHKERR mergeNodes(father, mother, bit, &level_tets, only_if_improve_quality,
484  move, th);
486 }
487 
488 } // namespace MoFEM
MoFEM::NodeMergerInterface::errorIfNoCommonEdge
bool errorIfNoCommonEdge
Send error if no common edge.
Definition: NodeMerger.hpp:128
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MortarCtestFunctions::check_tests
MoFEMErrorCode check_tests(int ss, int test_num, double &expected_energy, double &expected_contact_area, int &expected_nb_gauss_pts)
Definition: MortarCtestFunctions.hpp:229
MoFEM::NodeMergerInterface::query_interface
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Definition: NodeMerger.cpp:10
MoFEM::NodeMergerInterface::FaceMap
Definition: NodeMerger.hpp:160
MoFEM::NodeMergerInterface::FaceMapIdx
multi_index_container< FaceMap, indexed_by< hashed_unique< composite_key< FaceMap, member< FaceMap, EntityHandle, &FaceMap::n0 >, member< FaceMap, EntityHandle, &FaceMap::n1 > > > > > FaceMapIdx
Definition: NodeMerger.hpp:170
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
MoFEM::NodeMergerInterface::mergeNodes
MoFEMErrorCode mergeNodes(EntityHandle father, EntityHandle mother, Range &out_tets, Range *tets_ptr=NULL, const bool only_if_improve_quality=false, const double move=0, const int line_search=0, Tag th=NULL, const int verb=0)
merge nodes which sharing edge
Definition: NodeMerger.cpp:37
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::min_non_abs
static auto min_non_abs(const double a, const double b)
Definition: NodeMerger.cpp:16
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
MoFEM::NodeMergerInterface::NodeMergerInterface
NodeMergerInterface(const MoFEM::Core &core)
Definition: NodeMerger.cpp:20
MoFEM::NodeMergerInterface
Merge node by collapsing edge between them.
Definition: NodeMerger.hpp:18
THROW_MESSAGE
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:574
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::NodeMergerInterface::ParentChild
Definition: NodeMerger.hpp:102
MoFEM::NodeMergerInterface::minQualityFunction
boost::function< double(const double a, const double b)> minQualityFunction
Definition: NodeMerger.hpp:125
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::NodeMergerInterface::parentChildMap
ParentChildMap parentChildMap
Definition: NodeMerger.hpp:158
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
a
constexpr double a
Definition: approx_sphere.cpp:30
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
MoFEM::CoreTmp< 0 >::get_th_RefParentHandle
Tag get_th_RefParentHandle() const
Definition: Core.hpp:197
MoFEM::NodeMergerInterface::cOre
MoFEM::Core & cOre
Definition: NodeMerger.hpp:124
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::UnknownInterface
base class for all interface classes
Definition: UnknownInterface.hpp:34
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::NodeMergerInterface::lineSearch
MoFEMErrorCode lineSearch(Range &check_tests, EntityHandle father, EntityHandle mother, int line_search, FTensor::Tensor1< double, 3 > &t_move, Tag th=NULL)
Use bisection method to find point of edge collapse.
Definition: NodeMerger.cpp:403
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
UBlasVector< double >
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
MoFEM::NodeMergerInterface::getSubInterfaceOptions
MoFEMErrorCode getSubInterfaceOptions()
Definition: NodeMerger.cpp:24
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::Tools::volumeLengthQuality
static double volumeLengthQuality(const double *coords)
Calculate tetrahedron volume length quality.
Definition: Tools.cpp:15
MoFEM::NodeMergerInterface::successMerge
bool successMerge
True if marge is success.
Definition: NodeMerger.hpp:127
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
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
VERY_VERBOSE
@ VERY_VERBOSE
Definition: definitions.h:223
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MoFEM::BcManagerImplTools::get_adj_ents
auto get_adj_ents(moab::Interface &moab, const Range &ents)
Definition: BcManager.cpp:17
MoFEM::NodeMergerInterface::minQuality
MoFEMErrorCode minQuality(Range &check_tests, EntityHandle father, EntityHandle mother, double *coords_move, double &min_quality, Tag th=NULL, boost::function< double(double, double)> f=[](double a, double b) -> double { return std::min(a, b);})
Calualte quality if nodes merged.
Definition: NodeMerger.cpp:339