1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
|
#include "ff/Face.h"
#include "ff/Polyhedron.h"
#include "ff/Topology.h"
#include "test/3rdparty/catch.hpp"
#include <iostream>
#define CHECK_NEAR(a, b, eps) CHECK(fabs((a)-(b)) < eps);
// ************************************************************************************************
// Cartesian Tetrahedron
// ************************************************************************************************
ff::Polyhedron CartesianTetrahedron()
{
const ff::Topology topology = {
{{{0, 2, 1}, false}, {{0, 1, 3}, false}, {{0, 3, 2}, false}, {{1, 2, 3}, false}}, false};
constexpr double s = 1. / 4;
std::vector<R3> vertices{{-s, -s, -s}, {1 - s, -s, -s}, {-s, 1 - s, -s}, {-s, -s, 1 - s}};
return ff::Polyhedron(topology, vertices);
};
TEST_CASE("CartesianTetrahedron:Volume", "")
{
ff::Polyhedron P{CartesianTetrahedron()};
double volume = P.volume();
double volume2 = 0;
for (const ff::Face& face : P.faces())
volume2 += face.pyramidalVolume();
CHECK_NEAR(volume2, volume, 1E-13 * volume);
}
// ************************************************************************************************
// Pyramid4
// ************************************************************************************************
const ff::Topology topology = {{{{3, 2, 1, 0}, true},
{{0, 1, 5, 4}, false},
{{1, 2, 6, 5}, false},
{{2, 3, 7, 6}, false},
{{3, 0, 4, 7}, false},
{{4, 5, 6, 7}, true}},
false};
double r = .999;
double a = 1.;
double b = a * (1 - r);
double h = 3;
// center of mass: distance from bottom
double zcom = h * (1. / 4) * (6 - 8 * r + 3 * r * r) / (3 - 3 * r + r * r);
std::vector<R3> vertices{// base:
{-a, -a, -zcom},
{a, -a, -zcom},
{a, a, -zcom},
{-a, a, -zcom},
// top:
{-b, -b, h - zcom},
{b, -b, h - zcom},
{b, b, h - zcom},
{-b, b, h - zcom}};
ff::Polyhedron pyramid(topology, vertices);
TEST_CASE("Pyramid4:Volume", "")
{
double volume = pyramid.volume();
double volume2 = 0;
for (const ff::Face& face : pyramid.faces())
volume2 += face.pyramidalVolume();
CHECK_NEAR(volume2, volume, 1E-13 * volume);
}
TEST_CASE("Pyramid4:Bottom", "")
{
CHECK_NEAR(pyramid.z_bottom(), -zcom, 1E-13);
}
|