Father is sties, mother is merged.
Move node on the edge, 0 not move, 1 move to mother side, 0.5 will be in the middle.
58 CHKERR m_field.get_moab().get_adjacencies(&father, 1, 1,
false, father_edges);
60 CHKERR m_field.get_moab().get_adjacencies(&mother, 1, 1,
false, mother_edges);
64 CHKERR m_field.get_moab().get_adjacencies(&father, 1, 3,
false, father_tets);
66 CHKERR m_field.get_moab().get_adjacencies(&mother, 1, 3,
false, mother_tets);
67 if (tets_ptr != NULL) {
68 mother_tets = intersect(mother_tets, *tets_ptr);
69 father_tets = intersect(father_tets, *tets_ptr);
74 common_edge = intersect(father_edges, mother_edges);
75 if (tets_ptr != NULL) {
76 Range tets = unite(father_tets, mother_tets);
78 CHKERR m_field.get_moab().get_adjacencies(tets, 1,
false, tets_edges,
79 moab::Interface::UNION);
80 common_edge = intersect(common_edge, tets_edges);
81 father_edges = intersect(father_edges, tets_edges);
82 mother_edges = intersect(mother_edges, tets_edges);
88 "no common edge between nodes");
89 }
else if (common_edge.empty()) {
91 if (tets_ptr != NULL) {
92 seed_tets.merge(*tets_ptr);
101 CHKERR m_field.get_moab().get_adjacencies(common_edge, 3,
true, edge_tets);
103 if (tets_ptr != NULL) {
104 edge_tets = intersect(edge_tets, *tets_ptr);
107 mother_tets = subtract(mother_tets, edge_tets);
108 father_tets = subtract(father_tets, edge_tets);
111 const int num_nodes) {
114 CHKERR m_field.get_moab().get_coords(conn, num_nodes, &coords[0]);
116 CHKERR m_field.get_moab().tag_get_data(
th, conn, num_nodes, &coords[0]);
121 auto get_tensor = [](
VectorDouble &coords,
const int shift) {
123 &coords[shift], &coords[shift + 1], &coords[shift + 2]);
131 auto t_n0 = get_tensor(coords, 0);
132 auto t_n1 = get_tensor(coords, 3);
133 t_move(
i) = t_n0(
i) + move * (t_n1(
i) - t_n0(
i));
136 if (line_search > 0) {
137 Range
check_tests = unite(father_tets, mother_tets);
141 if (only_if_improve_quality) {
144 if (move > 0 || line_search) {
148 double min_quality0 = 1;
153 double min_quality = 1;
155 ((move > 0) || line_search) ? &t_move(0) : NULL,
157 if (min_quality < min_quality0) {
158 if (tets_ptr != NULL) {
159 out_tets = *tets_ptr;
167 if (move > 0 || line_search) {
169 CHKERR m_field.get_moab().set_coords(&father, 1, &t_move(0));
171 CHKERR m_field.get_moab().tag_set_data(
th, &father, 1, &t_move(0));
176 int *ret_num_nodes =
nullptr) {
179 CHKERR m_field.get_moab().get_connectivity(ent, conn, num_nodes,
true);
181 *ret_num_nodes = num_nodes;
185 auto create_tet = [
this, &m_field](
const EntityHandle *new_conn,
189 CHKERR m_field.get_moab().get_adjacencies(new_conn, 4, 3,
false, tets);
190 bool tet_found =
false;
191 for (
auto it_tet : tets) {
194 CHKERR m_field.get_moab().get_connectivity(it_tet, tet_conn, num_nodes,
196 const EntityHandle *
p = std::find(tet_conn, &tet_conn[4], new_conn[0]);
197 if (
p != tet_conn + 4) {
198 int s = std::distance(tet_conn,
p);
200 for (;
n != 4; ++
n) {
201 const int idx[] = {0, 1, 2, 3, 0, 1, 2, 3};
202 if (tet_conn[idx[s +
n]] != new_conn[
n])
205 if (
n == 4 && !tet_found) {
209 THROW_MESSAGE(
"More that one tet with the same connectivity");
215 CHKERR m_field.get_moab().create_element(MBTET, new_conn, 4, tet);
227 for (
auto m_tet : mother_tets) {
231 int nb_mother_verts = 0;
232 for (
int nn = 0; nn != 4; ++nn) {
233 if (conn[nn] == father) {
235 "Tet has father vertex, impossible but here it is");
237 if (conn[nn] == mother) {
238 new_conn[nn] = father;
241 new_conn[nn] = conn[nn];
244 if (nb_mother_verts != 1) {
246 "Tet should have only one vertex but have %d", nb_mother_verts);
252 created_tets.insert(create_tet(new_conn, m_tet));
256 auto get_adj_ents = [&](
const Range &ents,
const bool create) {
258 for (
int dd = 1;
dd <= 2;
dd++)
259 CHKERR m_field.get_moab().get_adjacencies(ents,
dd, create, adj,
260 moab::Interface::UNION);
263 auto adj_crated_ents = get_adj_ents(created_tets,
true);
264 adj_crated_ents.erase(common_edge[0]);
267 for (
auto ent : adj_crated_ents) {
272 bool father_node =
false;
273 for (
int nn = 0; nn != num_nodes; nn++) {
274 if (conn[nn] == father)
277 small_conn[ii++] = conn[nn];
281 std::sort(&small_conn[0], &small_conn[ii]);
283 face_map.insert(FaceMap(ent, small_conn[0], small_conn[1]));
285 face_map.insert(FaceMap(ent, small_conn[0], 0));
289 auto adj_mother_ents = get_adj_ents(mother_tets,
false);
290 adj_mother_ents.erase(common_edge[0]);
291 for (
auto ent : adj_mother_ents) {
298 for (; nn != num_nodes; ++nn) {
299 if (conn[nn] == mother) {
302 small_conn[ii++] = conn[nn];
305 if (nb_new_node > 0) {
307 std::sort(&small_conn[0], &small_conn[ii]);
313 FaceMapIdx::iterator fit = face_map.find(boost::make_tuple(n0, n1));
314 if (fit != face_map.end()) {
317 if (m_field.get_moab().dimension_from_handle(parent) !=
318 m_field.get_moab().dimension_from_handle(child))
320 "Huston we have a problem!");
333 if (tets_ptr != NULL) {
334 seed_tets.merge(*tets_ptr);
335 seed_tets = subtract(seed_tets, edge_tets);
336 seed_tets = subtract(seed_tets, mother_tets);
338 seed_tets.merge(created_tets);
339 out_tets.swap(seed_tets);
344 std::cerr <<
"Nodes merged" << endl;
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ MOFEM_DATA_INCONSISTENCY
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
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)
UBlasVector< double > VectorDouble
Tag get_th_RefParentHandle() const
multi_index_container< FaceMap, indexed_by< hashed_unique< composite_key< FaceMap, member< FaceMap, EntityHandle, &FaceMap::n0 >, member< FaceMap, EntityHandle, &FaceMap::n1 > > > > > FaceMapIdx
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.