1 #ifndef MEDUSA_BITS_DOMAINS_POLYHEDRONSHAPE_HPP_
2 #define MEDUSA_BITS_DOMAINS_POLYHEDRONSHAPE_HPP_
10 template <
typename vec_t>
12 const std::function<
scalar_t(vec_t)>& dr,
int type)
const {
13 bool vexists, fexists, cexists;
14 Mesh::Property_map<VI, Vector> vnormals;
15 Mesh::Property_map<FI, Vector> fnormals;
16 Mesh::Property_map<FI, CGAL::Color> fcolors;
17 std::tie(vnormals, vexists) = surface.property_map<
VI,
Vector>(
"v:normals");
18 assert_msg(vexists,
"Vertex normals were not computed.");
19 std::tie(fnormals, fexists) = surface.property_map<
FI,
Vector>(
"f:normals");
20 assert_msg(fexists,
"Face normals were not computed.");
21 std::tie(fcolors, cexists) = surface.property_map<
FI, CGAL::Color>(
"f:color");
22 auto get_face_type = [&](
FI face) {
23 return (cexists && type == 0) ?
24 rgb2type(fcolors[face].r(), fcolors[face].g(), fcolors[face].b()) :
29 auto add_point = [&](
const vec_t& p,
int type,
const vec_t& n) {
32 d.addBoundaryNode(p, type, n);
36 for (
VI idx : surface.vertices()) {
37 int face_type = get_face_type(surface.face(surface.halfedge(idx)));
38 auto n = vnormals[idx];
39 add_point(getPoint(idx), face_type, {n.x(), n.y(), n.z()});
42 for (Mesh::edge_index idx : surface.edges()) {
43 FI face1 = surface.face(idx.halfedge());
44 FI face2 = surface.face(surface.opposite(idx.halfedge()));
45 Vector normal1 = fnormals[face1];
46 Vector normal2 = fnormals[face2];
47 vec_t n1 = {normal1.x(), normal1.y(), normal1.z()};
48 vec_t n2 = {normal2.x(), normal2.y(), normal2.z()};
49 vec_t n = (n1 + n2).normalized();
50 int face_type = get_face_type(face1);
52 vec_t p1 = getPoint(surface.vertex(idx, 0));
53 vec_t p2 = getPoint(surface.vertex(idx, 1));
55 add_point(p, face_type, n);
59 for (
FI face : surface.faces()) {
60 auto range = surface.vertices_around_face(surface.halfedge(face));
61 assert_msg(range.size() == 3,
"All faces must be triangles.");
62 auto it = range.begin();
69 int face_type = get_face_type(face);
70 Vector normal = fnormals[face];
71 vec_t n = {normal.x(), normal.y(), normal.z()};
72 std::vector<vec_t> points_on_surface =
74 getPoint(p1), getPoint(p2), getPoint(p3), n, dr);
75 for (
const auto& p : points_on_surface) {
76 add_point(p, face_type, n);
84 #endif // MEDUSA_BITS_DOMAINS_POLYHEDRONSHAPE_HPP_