ExploringPhysicsWithCpp/Code/Chapter_13_rigid_body/main_chapter_13.cpp at master · zvanjak/ExploringPhysicsWithCpp · GitHub

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
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
#ifdef MML_USE_SINGLE_HEADER
#include "MML.h"
#else
#include "MMLBase.h"

#include "base/BaseUtils.h"
#include "base/VectorN.h"
#include "base/Geometry3D.h"

#include "core/CoordTransf.h"
#include "core/CoordTransf/CoordTransf3D.h"

#include "RigidBody/MomentOfInertiaCalculator.h"
#endif

using namespace MML;
using namespace MPL;

void Example13_tensor_of_inertia_discrete_masses()
{
	// define set of discrete masses
	double a = 1;
	Vector3Cartesian pos1(a, a, 0);
	Vector3Cartesian pos2(-a, a, 0);
	Vector3Cartesian pos3(-a, -a, 0);
	Vector3Cartesian pos4(a, -a, 0);
	Vector3Cartesian pos5(0, 0, 4 * a);
	double m1 = 2;
	double m2 = 1;
	double m3 = 4;
	double m4 = 1;
	double m5 = 1;

	DiscreteMass mass1(pos1, m1);
	DiscreteMass mass2(pos2, m2);
	DiscreteMass mass3(pos3, m3);
	DiscreteMass mass4(pos4, m4);
	DiscreteMass mass5(pos5, m5);

	std::vector<DiscreteMass> masses = { mass1, mass2, mass3, mass4, mass5 };
	DiscreteMassesConfig massesConfig(masses);

	DiscreteMassMomentOfInertiaTensorCalculator calculator(massesConfig);
	Tensor2<3> tensor_orig = calculator.calculate();

	std::cout << "Tensor of inertia: " << std::endl;
	std::cout << tensor_orig << std::endl;

	// investigating what happens if we change coord.system, in two cases:
	// 1. using coord.system transform we calculate TRANSFORMED (original) tensor
	// 2. using coord.system transform we calculate NEW set of masses and then calculate tensor

	// new coord.system is rotated around x axis for 30 degrees
	CoordTransfCart3DRotationXAxis coord_transf(Utils::DegToRad(30.0));

	// 1) - calculated tensor transformation
	Tensor2<3> tensor_transf = coord_transf.transfTensor2(tensor_orig, Vector3Cartesian(1, 1, 1));

	std::cout << "Tensor of inertia transformed: " << std::endl;
	std::cout << tensor_transf << std::endl;

	// 2) - change masses position and calculate new tensor
	DiscreteMassesConfig massesTransConfig(masses);
	for (auto& mass : massesTransConfig._masses)
		mass._position = coord_transf.transf(mass._position);

	DiscreteMassMomentOfInertiaTensorCalculator calculator2(massesTransConfig);
	Tensor2<3> tensor_changed = calculator2.calculate();

	std::cout << "Tensor of inertia rotated masses: " << std::endl;
	std::cout << tensor_changed << std::endl;
}

void Example13_tensor_of_inertia_continuous_mass()
{
	// create ContinuousMass representation for cube of constant density
	SolidBodyWithBoundaryConstDensity cube(-0.5, 0.5,
		[](Real x) { return -0.5; },
		[](Real x) { return 0.5; },
		[](Real x, Real y) { return -0.5; },
		[](Real x, Real y) { return 0.5; },
		1.0);

	ContinuousMassMomentOfInertiaTensorCalculator calculator(cube);

	Tensor2<3> tensor = calculator.calculate();

	std::cout << "Tensor of inertia for cube: " << std::endl;
	std::cout << tensor << std::endl;

	// let's do the same for sphere
	SolidBodyWithBoundaryConstDensity sphere(-1.0, 1.0,
		[](Real x) { return -sqrt(1 - x * x); },
		[](Real x) { return sqrt(1 - x * x); },
		[](Real x, Real y) { return -sqrt(1 - x * x - y * y); },
		[](Real x, Real y) { return sqrt(1 - x * x - y * y); },
		1.0);

	ContinuousMassMomentOfInertiaTensorCalculator calculator2(sphere);

	Tensor2<3> tensor2 = calculator2.calculate();

	std::cout << "Tensor of inertia for sphere: " << std::endl;
	std::cout << tensor2 << std::endl;

	std::cout << "Theoretical value for sphere: " << std::endl;
	std::cout << 2.0 / 5.0 * (4.0 / 3.0 * Constants::PI) << std::endl;
}

void Example13_tensor_of_inertia()
{
	std::cout << "***********************************************************************" << std::endl;
	std::cout << "****                   EXAMPLE 13 - tensor of inertia               ****" << std::endl;
	std::cout << "***********************************************************************" << std::endl;

	Example13_tensor_of_inertia_discrete_masses();
	Example13_tensor_of_inertia_continuous_mass();
}

int main()
{
	Example13_tensor_of_inertia();

  return 0;
}