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.
40 {
44
45
47 CHKERR m_field.get_moab().get_adjacencies(&father, 1, 1,
false, father_edges);
49 CHKERR m_field.get_moab().get_adjacencies(&mother, 1, 1,
false, mother_edges);
50
51
53 CHKERR m_field.get_moab().get_adjacencies(&father, 1, 3,
false, father_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
63 common_edge = intersect(father_edges, mother_edges);
64 if (tets_ptr != NULL) {
65 Range tets = unite(father_tets, mother_tets);
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
77 "no common edge between nodes");
78 } else if (common_edge.empty()) {
80 if (tets_ptr != NULL) {
81 seed_tets.merge(*tets_ptr);
82 }
83 out_tets = seed_tets;
86 }
87
88
90 CHKERR m_field.get_moab().get_adjacencies(common_edge, 3,
true, edge_tets);
91
92 if (tets_ptr != NULL) {
93 edge_tets = intersect(edge_tets, *tets_ptr);
94 }
95
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) {
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
117 if (move > 0) {
118 EntityHandle conn[] = {father, mother};
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) {
128 }
129
130 if (only_if_improve_quality) {
132
133 if (move > 0 || line_search) {
135 }
136
137 double min_quality0 = 1;
142 double min_quality = 1;
144 ((move > 0) || line_search) ? &t_move(0) : NULL,
146 if (min_quality < min_quality0) {
147 if (tets_ptr != NULL) {
148 out_tets = *tets_ptr;
149 }
152 }
153 }
154
155
156 if (move > 0 || line_search) {
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;
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);
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;
198 THROW_MESSAGE(
"More that one tet with the same connectivity");
199 }
200 }
201 }
202 if (!tet_found) {
203
204 CHKERR m_field.get_moab().create_element(MBTET, new_conn, 4, tet);
206 &tet, 1, &parent);
208 }
209 return tet;
210 };
211
212
214
216 for (auto m_tet : mother_tets) {
217 const EntityHandle *conn = get_conn(m_tet);
218 EntityHandle new_conn[4];
219
220 int nb_mother_verts = 0;
221 for (int nn = 0; nn != 4; ++nn) {
222 if (conn[nn] == father) {
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) {
235 "Tet should have only one vertex but have %d", nb_mother_verts);
236 }
237
239
240
241 created_tets.insert(create_tet(new_conn, m_tet));
242 }
243
244
245 auto get_adj_ents = [&](
const Range &ents,
const bool create) {
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
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))
309 "Huston we have a problem!");
310
311
313 &child, 1, &parent);
314
316 }
317 }
318 }
319
320
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
331
333 std::cerr << "Nodes merged" << endl;
334
336}
#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.
FTensor::Index< 'n', SPACE_DIM > n
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.
Public Member Functions inherited from MoFEM::UnknownInterface