@@ -44,6 +44,130 @@ void compute_curvature_measures(const data::Space<StructureT>& space, std::vecto
4444
4545 workspace.face_normals .resize (num_faces);
4646
47+ if constexpr (std::is_same_v<StructureT, data::DiscreteExteriorCalculus>) {
48+ const auto & x = space.x ;
49+ const auto & y = space.y ;
50+ const auto & z = space.z ;
51+ const auto & face_v0 = structure.face_v0 ;
52+ const auto & face_v1 = structure.face_v1 ;
53+ const auto & face_v2 = structure.face_v2 ;
54+ const auto & vertex_face_offsets = structure.vertex_face_offsets ;
55+ const auto & vertex_face_data = structure.vertex_face_data ;
56+
57+ core::parallel_for_index (
58+ 0 , static_cast <int >(num_faces),
59+ [&](int face_idx) {
60+ const size_t f = static_cast <size_t >(face_idx);
61+ const uint32_t i0 = face_v0[f];
62+ const uint32_t i1 = face_v1[f];
63+ const uint32_t i2 = face_v2[f];
64+
65+ const float e10x = x[i1] - x[i0];
66+ const float e10y = y[i1] - y[i0];
67+ const float e10z = z[i1] - z[i0];
68+ const float e20x = x[i2] - x[i0];
69+ const float e20y = y[i2] - y[i0];
70+ const float e20z = z[i2] - z[i0];
71+
72+ workspace.face_normals [f] = {
73+ e10x * e20y - e10y * e20x,
74+ e10y * e20z - e10z * e20y,
75+ e10z * e20x - e10x * e20z,
76+ };
77+ },
78+ 256 );
79+
80+ core::parallel_for_index (
81+ 0 , static_cast <int >(num_verts),
82+ [&](int vertex_idx) {
83+ const size_t i = static_cast <size_t >(vertex_idx);
84+ const uint32_t begin = vertex_face_offsets[i];
85+ const uint32_t end = vertex_face_offsets[i + 1 ];
86+ if (begin == end) {
87+ return ;
88+ }
89+
90+ float angle_sum = 0 .0f ;
91+ float area_sum = 0 .0f ;
92+
93+ float n_xy = 0 .0f ;
94+ float n_yz = 0 .0f ;
95+ float n_zx = 0 .0f ;
96+
97+ float sum_x = 0 .0f ;
98+ float sum_y = 0 .0f ;
99+ float sum_z = 0 .0f ;
100+
101+ const float px = x[i];
102+ const float py = y[i];
103+ const float pz = z[i];
104+
105+ for (uint32_t idx = begin; idx < end; ++idx) {
106+ const uint32_t f_idx = vertex_face_data[idx];
107+ const core::Bivec3 fn = workspace.face_normals [f_idx];
108+ n_xy += fn.xy ;
109+ n_yz += fn.yz ;
110+ n_zx += fn.zx ;
111+
112+ const uint32_t i0 = face_v0[f_idx];
113+ const uint32_t i1 = face_v1[f_idx];
114+ const uint32_t i2 = face_v2[f_idx];
115+
116+ uint32_t a = i0;
117+ uint32_t b = i1;
118+ if (i0 == i) {
119+ a = i1;
120+ b = i2;
121+ } else if (i1 == i) {
122+ a = i2;
123+ b = i0;
124+ }
125+
126+ const float ux = x[a] - px;
127+ const float uy = y[a] - py;
128+ const float uz = z[a] - pz;
129+ const float vx = x[b] - px;
130+ const float vy = y[b] - py;
131+ const float vz = z[b] - pz;
132+
133+ const float dot = ux * vx + uy * vy + uz * vz;
134+ const float wedge_xy = ux * vy - uy * vx;
135+ const float wedge_yz = uy * vz - uz * vy;
136+ const float wedge_zx = uz * vx - ux * vz;
137+ const float wedge_mag =
138+ std::sqrt (wedge_xy * wedge_xy + wedge_yz * wedge_yz + wedge_zx * wedge_zx);
139+
140+ angle_sum += std::atan2 (wedge_mag, dot);
141+ area_sum += 0 .5f * wedge_mag;
142+
143+ sum_x += x[a] + x[b];
144+ sum_y += y[a] + y[b];
145+ sum_z += z[a] + z[b];
146+ }
147+
148+ if (area_sum > 1e-12f ) {
149+ K[i] = static_cast <float >((2.0 * std::numbers::pi_v<double > - angle_sum) /
150+ (static_cast <double >(area_sum) / 3.0 ));
151+ }
152+
153+ const float n_mag_sq = n_xy * n_xy + n_yz * n_yz + n_zx * n_zx;
154+ const float n_inv = (n_mag_sq > 1e-12f ) ? 1 .0f / std::sqrt (n_mag_sq) : 0 .0f ;
155+
156+ const float normal_x = n_yz * n_inv;
157+ const float normal_y = n_zx * n_inv;
158+ const float normal_z = n_xy * n_inv;
159+
160+ const float inv_c = 0 .5f / static_cast <float >(end - begin);
161+ const float laplacian_x = sum_x * inv_c - px;
162+ const float laplacian_y = sum_y * inv_c - py;
163+ const float laplacian_z = sum_z * inv_c - pz;
164+
165+ H[i] = laplacian_x * normal_x + laplacian_y * normal_y + laplacian_z * normal_z;
166+ },
167+ 128 );
168+ return ;
169+ }
170+
47171 const auto get_face_vertex = [&](size_t face_idx, int corner) -> uint32_t {
48172 if constexpr (std::is_same_v<StructureT, data::DiscreteExteriorCalculus>) {
49173 if (corner == 0 )
0 commit comments