12 {
14
15 Range fineFeatureTris;
19 if (vols.empty()) {
21 "no 3d entities found in the mesh");
22 }
23
24 Skinner skinner(&m_field.
get_moab());
26 CHKERR skinner.find_skin(0, vols,
false, skin_tris);
27 if (skin_tris.empty()) {
29 "no skin triangles found in the mesh");
30 }
31
36 normal.normalize();
38 };
39
40 const double fine_planar_cos = std::cos(1e-6);
41
43 CHKERR m_field.
get_moab().get_adjacencies(skin_tris, 1,
true, skin_edges,
44 moab::Interface::UNION);
45
46 for (auto edge : skin_edges) {
48 CHKERR m_field.
get_moab().get_adjacencies(&edge, 1, 2,
false, adj_tris,
49 moab::Interface::UNION);
50 adj_tris = intersect(adj_tris, skin_tris);
51 if (adj_tris.size() != 2) {
52 continue;
53 }
54
57 CHKERR tri_normal(adj_tris[0], n0);
58 CHKERR tri_normal(adj_tris[1], n1);
59 const double dot = n0(
i) * n1(
i);
60 if (std::abs(dot) <= fine_planar_cos) {
61 fineFeatureTris.insert(adj_tris[0]);
62 fineFeatureTris.insert(adj_tris[1]);
63 output_edges.insert(edge);
64 }
65 }
66
67 std::map<EntityHandle, std::vector<EntityHandle>> vert_edges;
68 for (auto edge : output_edges) {
70 int nconn = 0;
71 CHKERR m_field.
get_moab().get_connectivity(edge, conn, nconn,
false);
72 if (nconn != 2)
73 continue;
74 vert_edges[conn[0]].push_back(edge);
75 vert_edges[conn[1]].push_back(edge);
76 }
77
79 for (auto edge : output_edges) {
81 int nconn = 0;
82 CHKERR m_field.
get_moab().get_connectivity(edge, conn, nconn,
false);
83 if (nconn != 2)
84 continue;
85
86 bool keep_edge = true;
87 for (int vv = 0; vv < 2; ++vv) {
88 const auto v = conn[vv];
89 const auto it = vert_edges.find(
v);
90 if (it == vert_edges.end()) {
91 keep_edge = false;
92 break;
93 }
94 const auto &inc_edges = it->second;
95 if (inc_edges.size() > 2) {
96 keep_edge = false;
97 break;
98 }
99 if (inc_edges.size() == 2) {
100 const auto other_edge =
101 (inc_edges[0] == edge) ? inc_edges[1] : inc_edges[0];
103 int nconn_other = 0;
104 CHKERR m_field.
get_moab().get_connectivity(other_edge, conn_other,
105 nconn_other, false);
106 if (nconn_other != 2) {
107 keep_edge = false;
108 break;
109 }
111 (conn_other[0] ==
v) ? conn_other[1] : conn_other[0];
112
113 double pv[3], pa[3], pb[3];
118
122 &pv[0], &pv[1], &pv[2]);
124 &pa[0], &pa[1], &pa[2]);
126 &pb[0], &pb[1], &pb[2]);
127 a(
i) = t_pa(
i) - t_pv(
i);
128 b(
i) = t_pb(
i) - t_pv(
i);
129 const double na =
a.l2();
130 const double nb = b.
l2();
131 if (na == 0 || nb == 0) {
132 keep_edge = false;
133 break;
134 }
135 const double dot = (
a(
i) * b(
i)) / (na * nb);
136 if (std::abs(dot) < fine_planar_cos) {
137 keep_edge = false;
138 break;
139 }
140 }
141 }
142
143 if (keep_edge) {
144 straight_edges.insert(edge);
145 }
146 }
147
148 output_edges = straight_edges;
149 output_vertices.clear();
150 if (!output_edges.empty()) {
151 std::map<EntityHandle, int> vert_count;
152 for (auto edge : output_edges) {
154 int nconn = 0;
155 CHKERR m_field.
get_moab().get_connectivity(edge, conn, nconn,
false);
156 if (nconn != 2)
157 continue;
158 vert_count[conn[0]]++;
159 vert_count[conn[1]]++;
160 }
161 Range filtered_edges;
162 for (auto edge : output_edges) {
164 int nconn = 0;
165 CHKERR m_field.
get_moab().get_connectivity(edge, conn, nconn,
false);
166 if (nconn != 2)
167 continue;
168 const int d0 = vert_count[conn[0]];
169 const int d1 = vert_count[conn[1]];
170 if (!(d0 == 1 && d1 == 1)) {
171 filtered_edges.insert(edge);
172 for (int vv = 0; vv < 2; ++vv) {
173 const auto v = conn[vv];
174 auto it = vert_count.find(
v);
175 if (it != vert_count.end() && it->second == 1) {
176 output_vertices.insert(
v);
177 }
178 }
179 }
180 }
181 output_edges = filtered_edges;
182 }
183
184 if (!fineFeatureTris.empty())
186 "tri_features_test.vtk", fineFeatureTris);
187 if (!output_edges.empty())
189 "edges_test.vtk", output_edges);
190 if (!output_vertices.empty())
192 "vertices_test.vtk", output_vertices);
193
195 }
#define FTENSOR_INDEX(DIM, I)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
virtual moab::Interface & get_moab()=0
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
