Wildmeshing Toolkit
Loading...
Searching...
No Matches
marching.cpp
Go to the documentation of this file.
1#include <catch2/catch_test_macros.hpp>
2
3#include <tools/DEBUG_TetMesh.hpp>
4#include <tools/DEBUG_TriMesh.hpp>
5#include <tools/TetMesh_examples.hpp>
6#include <tools/TriMesh_examples.hpp>
7
8#include <wmtk/io/MeshReader.hpp>
9#include <wmtk/io/ParaviewWriter.hpp>
10#include <wmtk/simplex/link.hpp>
11#include <wmtk/utils/mesh_utils.hpp>
12
13#include <wmtk/components/marching/marching.hpp>
14
15using namespace wmtk;
16
17
18TEST_CASE("marching_component_tri", "[components][marching]")
19{
20 const int64_t input_tag_value_0 = 0;
21 const int64_t input_tag_value_1 = 1;
22 const int64_t isosurface_tag_value = 2;
23
24 // 0---1---2
25 // / \ / \ / \ .
26 // 3---4---5---6
27 // \ / \ / .
28 // 7---8
29 tests::DEBUG_TriMesh m = tests::hex_plus_two_with_position();
30
31 components::MarchingOptions options;
32 options.position_handle = m.get_attribute_handle<double>("vertices", PrimitiveType::Vertex);
33 options.label_handles[PrimitiveType::Vertex] = m.register_attribute<int64_t>(
34 "vertex_tag",
35 PrimitiveType::Vertex,
36 1,
37 false,
38 input_tag_value_0);
39 options.label_handles[PrimitiveType::Edge] =
40 m.register_attribute<int64_t>("marching_edge_tag", PrimitiveType::Edge, 1);
41
42 options.input_values = {input_tag_value_0, input_tag_value_1};
43 options.output_value = isosurface_tag_value;
44
45 int64_t expected_isosurface_vertex_num = 0;
46 int64_t expected_isosurface_edge_num = 0;
47
48 SECTION("4")
49 {
50 const std::vector<Tuple>& vertex_tuples = m.get_all(wmtk::PrimitiveType::Vertex);
51 attribute::Accessor<int64_t> acc_vertex_tag =
52 m.create_accessor<int64_t>(options.label_handles[PrimitiveType::Vertex]);
53 acc_vertex_tag.scalar_attribute(vertex_tuples[4]) = input_tag_value_1;
54
55 expected_isosurface_vertex_num = 6;
56 expected_isosurface_edge_num = 6;
57 }
58 SECTION("4-5")
59 {
60 const std::vector<Tuple>& vertex_tuples = m.get_all(wmtk::PrimitiveType::Vertex);
61 attribute::Accessor<int64_t> acc_vertex_tag =
62 m.create_accessor<int64_t>(options.label_handles[PrimitiveType::Vertex]);
63 acc_vertex_tag.scalar_attribute(vertex_tuples[4]) = input_tag_value_1;
64 acc_vertex_tag.scalar_attribute(vertex_tuples[5]) = input_tag_value_1;
65
66 expected_isosurface_vertex_num = 9;
67 expected_isosurface_edge_num = 8;
68 }
69 SECTION("0-4-5")
70 {
71 const std::vector<Tuple>& vertex_tuples = m.get_all(wmtk::PrimitiveType::Vertex);
72 attribute::Accessor<int64_t> acc_vertex_tag =
73 m.create_accessor<int64_t>(options.label_handles[PrimitiveType::Vertex]);
74 acc_vertex_tag.scalar_attribute(vertex_tuples[0]) = input_tag_value_1;
75 acc_vertex_tag.scalar_attribute(vertex_tuples[4]) = input_tag_value_1;
76 acc_vertex_tag.scalar_attribute(vertex_tuples[5]) = input_tag_value_1;
77
78 expected_isosurface_vertex_num = 10;
79 expected_isosurface_edge_num = 8;
80 }
81 SECTION("0-4-5-with-filter")
82 {
83 const std::vector<Tuple>& vertex_tuples = m.get_all(wmtk::PrimitiveType::Vertex);
84 attribute::Accessor<int64_t> acc_vertex_tag =
85 m.create_accessor<int64_t>(options.label_handles[PrimitiveType::Vertex]);
86 acc_vertex_tag.scalar_attribute(vertex_tuples[0]) = input_tag_value_1;
87 acc_vertex_tag.scalar_attribute(vertex_tuples[4]) = input_tag_value_1;
88 acc_vertex_tag.scalar_attribute(vertex_tuples[5]) = input_tag_value_1;
89
90 attribute::MeshAttributeHandle filter =
91 m.register_attribute<int64_t>("edge_filter", PrimitiveType::Edge, 1);
92 options.edge_filter_handles.emplace_back(filter, 1);
93
94 attribute::Accessor<int64_t> acc_filter = m.create_accessor<int64_t>(filter);
95 acc_filter.scalar_attribute(m.edge_tuple_from_vids(0, 1)) = 1;
96 acc_filter.scalar_attribute(m.edge_tuple_from_vids(1, 4)) = 1;
97 acc_filter.scalar_attribute(m.edge_tuple_from_vids(1, 5)) = 1;
98 acc_filter.scalar_attribute(m.edge_tuple_from_vids(2, 5)) = 1;
99 acc_filter.scalar_attribute(m.edge_tuple_from_vids(5, 6)) = 1;
100
101 expected_isosurface_vertex_num = 5;
102 expected_isosurface_edge_num = 4;
103 }
104
105 int64_t expected_vertex_num =
106 m.get_all(PrimitiveType::Vertex).size() + expected_isosurface_vertex_num;
107
108 components::marching(m, options);
109
110 const auto& vertices = m.get_all(PrimitiveType::Vertex);
111 attribute::Accessor<int64_t> acc_vertex_tag =
112 m.create_accessor<int64_t>(options.label_handles[PrimitiveType::Vertex]);
113 // vertex number should be correct
114 {
115 CHECK(vertices.size() == expected_vertex_num);
116
117 int64_t isosurface_vertex_num = 0;
118 for (const Tuple& v : vertices) {
119 if (acc_vertex_tag.scalar_attribute(v) == isosurface_tag_value) {
120 isosurface_vertex_num++;
121 }
122 }
123 CHECK(isosurface_vertex_num == expected_isosurface_vertex_num);
124 }
125
126 const auto& edges = m.get_all(PrimitiveType::Edge);
127 wmtk::attribute::Accessor<int64_t> acc_edge_tag =
128 m.create_accessor<int64_t>(options.label_handles[PrimitiveType::Edge]);
129 // edge number should be correct
130 {
131 int64_t isosurface_edge_num = 0;
132 for (const Tuple& e : edges) {
133 if (acc_edge_tag.scalar_attribute(e) == isosurface_tag_value) {
134 isosurface_edge_num++;
135 }
136 }
137 CHECK(isosurface_edge_num == expected_isosurface_edge_num);
138 }
139
140 // should be manifold
141 {
142 for (const Tuple& v : vertices) {
143 if (acc_vertex_tag.scalar_attribute(v) == isosurface_tag_value) {
144 std::vector<Tuple> one_ring = simplex::link(m, simplex::Simplex::vertex(m, v))
145 .simplex_vector_tuples(PrimitiveType::Vertex);
146
147 int64_t tagged_neighbors = 0;
148 for (const Tuple& neigh : one_ring) {
149 if (acc_vertex_tag.scalar_attribute(neigh) == isosurface_tag_value) {
150 ++tagged_neighbors;
151 }
152 }
153
154 if (m.is_boundary_vertex(v)) {
155 CHECK(tagged_neighbors == 1);
156 } else {
157 CHECK(tagged_neighbors == 2);
158 }
159 }
160 }
161 }
162
163 if (false) {
164 wmtk::io::ParaviewWriter
165 writer("marching_2d_result", "vertices", m, true, true, true, false);
166 m.serialize(writer);
167 }
168}
169
170TEST_CASE("marching_component_tet", "[components][marching]")
171{
172 const int64_t input_tag_value_0 = 0;
173 const int64_t input_tag_value_1 = 1;
174 const int64_t isosurface_tag_value = 2;
175
176 // 0 ---------- 4
177 // / \\ // \ .
178 // / \ \ // \ .
179 // / \ \ // \ .
180 // / \ \3 \ .
181 // 1 --------- 2/ -------- 5 tuple edge 2-3
182 //
183
184 tests_3d::DEBUG_TetMesh m = tests_3d::three_incident_tets_with_positions();
185
186 components::MarchingOptions options;
187 options.position_handle = m.get_attribute_handle<double>("vertices", PrimitiveType::Vertex);
188 options.label_handles[PrimitiveType::Vertex] = m.register_attribute<int64_t>(
189 "vertex_tag",
190 PrimitiveType::Vertex,
191 1,
192 false,
193 input_tag_value_0);
194 options.label_handles[PrimitiveType::Edge] =
195 m.register_attribute<int64_t>("marching_edge_tag", PrimitiveType::Edge, 1);
196 options.label_handles[PrimitiveType::Triangle] =
197 m.register_attribute<int64_t>("marching_face_tag", PrimitiveType::Triangle, 1);
198
199 options.input_values = {input_tag_value_0, input_tag_value_1};
200 options.output_value = isosurface_tag_value;
201
202
203 int64_t expected_isosurface_vertex_num = 0;
204 int64_t expected_isosurface_face_num = 0;
205
206 SECTION("2-3")
207 {
208 const std::vector<Tuple>& vertex_tuples = m.get_all(wmtk::PrimitiveType::Vertex);
209 attribute::Accessor<int64_t> acc_vertex_tag =
210 m.create_accessor<int64_t>(options.label_handles[PrimitiveType::Vertex]);
211 acc_vertex_tag.scalar_attribute(vertex_tuples[2]) = input_tag_value_1;
212 acc_vertex_tag.scalar_attribute(vertex_tuples[3]) = input_tag_value_1;
213
214 expected_isosurface_vertex_num = 8;
215 expected_isosurface_face_num = 6;
216 }
217 SECTION("4-5")
218 {
219 const std::vector<Tuple>& vertex_tuples = m.get_all(wmtk::PrimitiveType::Vertex);
220 wmtk::attribute::Accessor<int64_t> acc_vertex_tag =
221 m.create_accessor<int64_t>(options.label_handles[PrimitiveType::Vertex]);
222 acc_vertex_tag.scalar_attribute(vertex_tuples[4]) = input_tag_value_1;
223 acc_vertex_tag.scalar_attribute(vertex_tuples[5]) = input_tag_value_1;
224
225 expected_isosurface_vertex_num = 5;
226 expected_isosurface_face_num = 3;
227 }
228 SECTION("1-0-4-5")
229 {
230 const std::vector<Tuple>& vertex_tuples = m.get_all(wmtk::PrimitiveType::Vertex);
231 wmtk::attribute::Accessor<int64_t> acc_vertex_tag =
232 m.create_accessor<int64_t>(options.label_handles[PrimitiveType::Vertex]);
233 acc_vertex_tag.scalar_attribute(vertex_tuples[1]) = input_tag_value_1;
234 acc_vertex_tag.scalar_attribute(vertex_tuples[0]) = input_tag_value_1;
235 acc_vertex_tag.scalar_attribute(vertex_tuples[4]) = input_tag_value_1;
236 acc_vertex_tag.scalar_attribute(vertex_tuples[5]) = input_tag_value_1;
237
238 expected_isosurface_vertex_num = 8;
239 expected_isosurface_face_num = 6;
240 }
241
242 int64_t expected_vertex_num =
243 m.get_all(PrimitiveType::Vertex).size() + expected_isosurface_vertex_num;
244
245
246 std::vector<attribute::MeshAttributeHandle> pass_through_attributes;
247 pass_through_attributes.emplace_back(
248 m.get_attribute_handle<double>("vertices", PrimitiveType::Vertex));
249
250 components::marching(m, options);
251
252 const auto& vertices = m.get_all(PrimitiveType::Vertex);
253 attribute::Accessor<int64_t> acc_vertex_tag =
254 m.create_accessor<int64_t>(options.label_handles[PrimitiveType::Vertex]);
255 // vertex number should be correct
256 {
257 CHECK(vertices.size() == expected_vertex_num);
258
259 int64_t isosurface_vertex_num = 0;
260 for (const Tuple& v : vertices) {
261 if (acc_vertex_tag.scalar_attribute(v) == isosurface_tag_value) {
262 isosurface_vertex_num++;
263 }
264 }
265 CHECK(isosurface_vertex_num == expected_isosurface_vertex_num);
266 }
267
268 // face number should be correct
269 const auto& faces = m.get_all(PrimitiveType::Triangle);
270 wmtk::attribute::Accessor<int64_t> acc_face_tag =
271 m.create_accessor<int64_t>(options.label_handles[PrimitiveType::Triangle]);
272 {
273 int64_t isosurface_face_num = 0;
274 for (const Tuple& f : faces) {
275 if (acc_face_tag.scalar_attribute(f) == isosurface_tag_value) {
276 isosurface_face_num++;
277 }
278 }
279 CHECK(isosurface_face_num == expected_isosurface_face_num);
280 }
281
282 if (false) {
283 wmtk::io::ParaviewWriter
284 writer("marching_3d_result", "vertices", m, true, true, true, true);
285 m.serialize(writer);
286 }
287}
wmtk::Tuple
The Tuple is the basic navigation tool in our mesh data structure.
Definition Tuple.hpp:19
wmtk::attribute::Accessor
A CachingAccessor that uses tuples for accessing attributes instead of indices.
Definition Accessor.hpp:28
wmtk::attribute::Accessor::scalar_attribute
T & scalar_attribute(const ArgType &t)
wmtk::simplex::SimplexCollection::simplex_vector_tuples
std::vector< Tuple > simplex_vector_tuples(PrimitiveType ptype) const
Return vector of all simplices of the requested type, as tuples.
Definition SimplexCollection.cpp:28
wmtk::simplex::Simplex::vertex
static Simplex vertex(const Mesh &m, const Tuple &t)
Definition Simplex.hpp:56
mesh_utils.hpp
wmtk::components::marching
void marching(Mesh &mesh, const MarchingOptions &options)
Perform maching tetrahedra/triangles.
Definition marching.cpp:11
wmtk::simplex::link
SimplexCollection link(const Mesh &mesh, const simplex::Simplex &simplex, const bool sort_and_clean)
Definition link.cpp:84
wmtk::components::MarchingOptions::edge_filter_handles
std::vector< std::pair< attribute::MeshAttributeHandle, int64_t > > edge_filter_handles
Filters (int64_t) on edge labels.
Definition MarchingOptions.hpp:25
wmtk::components::MarchingOptions::position_handle
attribute::MeshAttributeHandle position_handle
vertex positions (double)
Definition MarchingOptions.hpp:13
wmtk::components::MarchingOptions::input_values
std::vector< int64_t > input_values
The label values (int64_t) in between the marching surface should be placed.
Definition MarchingOptions.hpp:36
wmtk::components::MarchingOptions::output_value
int64_t output_value
The output label value (int 64_t) is assigned to the marching surface in all primitive types for whic...
Definition MarchingOptions.hpp:41
wmtk::components::MarchingOptions::label_handles
std::map< PrimitiveType, attribute::MeshAttributeHandle > label_handles
Labels for inside, outside, and marching surface (int64_t).
Definition MarchingOptions.hpp:20
TEST_CASE
TEST_CASE("marching_component_tri", "[components][marching]")
Definition marching.cpp:18