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.
37 {
41
42
44 CHKERR m_field.get_moab().get_adjacencies(&father, 1, 1,
false, father_edges);
46 CHKERR m_field.get_moab().get_adjacencies(&mother, 1, 1,
false, mother_edges);
47
48
50 CHKERR m_field.get_moab().get_adjacencies(&father, 1, 3,
false, father_tets);
52 CHKERR m_field.get_moab().get_adjacencies(&mother, 1, 3,
false, mother_tets);
53 if (tets_ptr != NULL) {
54 mother_tets = intersect(mother_tets, *tets_ptr);
55 father_tets = intersect(father_tets, *tets_ptr);
56 }
57
58
60 common_edge = intersect(father_edges, mother_edges);
61 if (tets_ptr != NULL) {
62 Range tets = unite(father_tets, mother_tets);
64 CHKERR m_field.get_moab().get_adjacencies(tets, 1,
false, tets_edges,
65 moab::Interface::UNION);
66 common_edge = intersect(common_edge, tets_edges);
67 father_edges = intersect(father_edges, tets_edges);
68 mother_edges = intersect(mother_edges, tets_edges);
69 }
70
71
74 "no common edge between nodes");
75 } else if (common_edge.empty()) {
77 if (tets_ptr != NULL) {
78 seed_tets.merge(*tets_ptr);
79 }
80 out_tets = seed_tets;
83 }
84
85
87 CHKERR m_field.get_moab().get_adjacencies(common_edge, 3,
true, edge_tets);
88
89 if (tets_ptr != NULL) {
90 edge_tets = intersect(edge_tets, *tets_ptr);
91 }
92
93 mother_tets = subtract(mother_tets, edge_tets);
94 father_tets = subtract(father_tets, edge_tets);
95
97 const int num_nodes) {
98 VectorDouble coords(3 * num_nodes);
100 CHKERR m_field.get_moab().get_coords(conn, num_nodes, &coords[0]);
101 } else {
102 CHKERR m_field.get_moab().tag_get_data(
th, conn, num_nodes, &coords[0]);
103 }
104 return coords;
105 };
106
107 auto get_tensor = [](VectorDouble &coords, const int shift) {
109 &coords[shift], &coords[shift + 1], &coords[shift + 2]);
110 };
111
112
114 if (move > 0) {
116 VectorDouble coords = get_coords(
th, conn, 2);
117 auto t_n0 = get_tensor(coords, 0);
118 auto t_n1 = get_tensor(coords, 3);
119 t_move(
i) = t_n0(
i) + move * (t_n1(
i) - t_n0(
i));
120 }
121
122 if (line_search > 0) {
123 Range check_tests = unite(father_tets, mother_tets);
125 }
126
127 if (only_if_improve_quality) {
128 Range check_tests = mother_tets;
129
130 if (move > 0 || line_search) {
131 check_tests.merge(father_tets);
132 }
133
134 double min_quality0 = 1;
139 double min_quality = 1;
141 ((move > 0) || line_search) ? &t_move(0) : NULL,
143 if (min_quality < min_quality0) {
144 if (tets_ptr != NULL) {
145 out_tets = *tets_ptr;
146 }
149 }
150 }
151
152
153 if (move > 0 || line_search) {
155 CHKERR m_field.get_moab().set_coords(&father, 1, &t_move(0));
156 } else {
157 CHKERR m_field.get_moab().tag_set_data(
th, &father, 1, &t_move(0));
158 }
159 }
160
162 int *ret_num_nodes = nullptr) {
163 int num_nodes;
165 CHKERR m_field.get_moab().get_connectivity(ent, conn, num_nodes,
true);
166 if (ret_num_nodes)
167 *ret_num_nodes = num_nodes;
168 return conn;
169 };
170
171 auto create_tet = [
this, &m_field](
const EntityHandle *new_conn,
175 CHKERR m_field.get_moab().get_adjacencies(new_conn, 4, 3,
false, tets);
176 bool tet_found = false;
177 for (auto it_tet : tets) {
179 int num_nodes;
180 CHKERR m_field.get_moab().get_connectivity(it_tet, tet_conn, num_nodes,
181 true);
182 const EntityHandle *p = std::find(tet_conn, &tet_conn[4], new_conn[0]);
183 if (p != tet_conn + 4) {
184 int s = std::distance(tet_conn, p);
186 for (;
n != 4; ++
n) {
187 const int idx[] = {0, 1, 2, 3, 0, 1, 2, 3};
188 if (tet_conn[idx[s +
n]] != new_conn[
n])
189 break;
190 }
191 if (
n == 4 && !tet_found) {
192 tet = it_tet;
193 tet_found = true;
195 THROW_MESSAGE(
"More that one tet with the same connectivity");
196 }
197 }
198 }
199 if (!tet_found) {
200
201 CHKERR m_field.get_moab().create_element(MBTET, new_conn, 4, tet);
203 &tet, 1, &parent);
205 }
206 return tet;
207 };
208
209
211
213 for (auto m_tet : mother_tets) {
216
217 int nb_mother_verts = 0;
218 for (int nn = 0; nn != 4; ++nn) {
219 if (conn[nn] == father) {
221 "Tet has father vertex, impossible but here it is");
222 }
223 if (conn[nn] == mother) {
224 new_conn[nn] = father;
225 nb_mother_verts++;
226 } else {
227 new_conn[nn] = conn[nn];
228 }
229 }
230 if (nb_mother_verts != 1) {
232 "Tet should have only one vertex but have %d", nb_mother_verts);
233 }
234
235 VectorDouble new_coords = get_coords(
th, new_conn, 4);
236
237
238 created_tets.insert(create_tet(new_conn, m_tet));
239 }
240
241
244 for (
int dd = 1;
dd <= 2;
dd++)
245 CHKERR m_field.get_moab().get_adjacencies(ents, dd, create, adj,
246 moab::Interface::UNION);
247 return adj;
248 };
249 auto adj_crated_ents =
get_adj_ents(created_tets,
true);
250 adj_crated_ents.erase(common_edge[0]);
251
253 for (auto ent : adj_crated_ents) {
254 int num_nodes;
257 int ii = 0;
258 bool father_node = false;
259 for (int nn = 0; nn != num_nodes; nn++) {
260 if (conn[nn] == father)
261 father_node = true;
262 else
263 small_conn[ii++] = conn[nn];
264 }
265 if (father_node) {
266 if (ii > 1)
267 std::sort(&small_conn[0], &small_conn[ii]);
268 if (ii == 2)
269 face_map.insert(FaceMap(ent, small_conn[0], small_conn[1]));
270 else
271 face_map.insert(FaceMap(ent, small_conn[0], 0));
272 }
273 }
274
275 auto adj_mother_ents =
get_adj_ents(mother_tets,
false);
276 adj_mother_ents.erase(common_edge[0]);
277 for (auto ent : adj_mother_ents) {
278 int num_nodes;
281 int nb_new_node = 0;
282 int nn = 0;
283 int ii = 0;
284 for (; nn != num_nodes; ++nn) {
285 if (conn[nn] == mother) {
286 nb_new_node++;
287 } else {
288 small_conn[ii++] = conn[nn];
289 }
290 }
291 if (nb_new_node > 0) {
292 if (ii > 1)
293 std::sort(&small_conn[0], &small_conn[ii]);
294
296 if (ii == 2)
297 n1 = small_conn[1];
298
299 FaceMapIdx::iterator fit = face_map.find(boost::make_tuple(n0, n1));
300 if (fit != face_map.end()) {
303 if (m_field.get_moab().dimension_from_handle(parent) !=
304 m_field.get_moab().dimension_from_handle(child))
306 "Huston we have a problem!");
307
308
310 &child, 1, &parent);
311
313 }
314 }
315 }
316
317
319 if (tets_ptr != NULL) {
320 seed_tets.merge(*tets_ptr);
321 seed_tets = subtract(seed_tets, edge_tets);
322 seed_tets = subtract(seed_tets, mother_tets);
323 }
324 seed_tets.merge(created_tets);
325 out_tets.swap(seed_tets);
326
328
330 std::cerr << "Nodes merged" << endl;
331
333}
#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 double n
refractive index of diffusive medium
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)
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.