GSI_Toolbox 1.0.0 A toolbox for Gas-Surface Interaction simulations

Loading...
Searching...
No Matches
Geometry.h
1
10template <typename T>
11 class Geometry {
12
13 private:
14
15 unsigned int num_vertices;
16 char header[81];
17 Matrix<T> p1[4], p2[4], p3[4];
18 Matrix<T> norm_x, norm_y, norm_z, norm_mod;
19 T A_ref, c_ref;
20 std::vector<Matrix<T>> GSI_prop;
21 std::vector<Matrix<T>> surface_prop;
22 std::vector<const char*> GSI_prop_names;
23 std::vector<const char*> surface_prop_names;
24 Matrix<T> bounding_box;
25
26 public:
27
33 Geometry();
34
40 Geometry(std::string filename);
41
47 ~Geometry();
48
54 unsigned int get_num_vertices() const { return num_vertices; }
55
61 void shift(T* offset);
62
68 void shift(Matrix<T> offset);
69
75 void scale(T factor);
76
83 void rotate(int axis, T angle);
84
91 void add_GSI_property(const char* property_name, T value);
92
99 void add_GSI_property(const char* property_name, const Matrix<T>& values);
100
107 void add_surface_property(const char* property_name, T value);
108
115 void add_surface_property(const char* property_name, const Matrix<T>& values);
116
122 void set_A_ref(T A_ref) { this->A_ref = A_ref; }
123
129 void set_c_ref(T c_ref) { this->c_ref = c_ref; }
130
134 void compute_bounding_box();
135
141 T get_A_ref() { return A_ref; }
142
148 T get_c_ref() { return c_ref; }
149
156 T compute_area(unsigned int index);
157
161 void print();
162
163 // Getter methods for the vertex coordinates and normal vectors.
164 Matrix<T>& get_p1_x() { return this->p1[0]; }
165 Matrix<T>& get_p1_y() { return this->p1[1]; }
166 Matrix<T>& get_p1_z() { return this->p1[2]; }
167 Matrix<T>& get_p2_x() { return this->p2[0]; }
168 Matrix<T>& get_p2_y() { return this->p2[1]; }
169 Matrix<T>& get_p2_z() { return this->p2[2]; }
170 Matrix<T>& get_p3_x() { return this->p3[0]; }
171 Matrix<T>& get_p3_y() { return this->p3[1]; }
172 Matrix<T>& get_p3_z() { return this->p3[2]; }
173 Matrix<T>& get_norm_x() { return this->norm_x; }
174 Matrix<T>& get_norm_y() { return this->norm_y; }
175 Matrix<T>& get_norm_z() { return this->norm_z; }
176
183 Matrix<T>& get_GSI_property(int index) { return GSI_prop[index]; }
184
191 char* get_GSI_property_name(int index) { return GSI_prop_names[index]; }
192
198 Matrix<T> get_bounding_box();
199
205 Matrix<T> get_center();
206
212 std::vector<Matrix<T>>& get_GSI_properties() { return GSI_prop; }
213
219 std::vector<const char*>& get_GSI_property_names() { return GSI_prop_names; }
220
226 std::vector<Matrix<T>>& get_surface_properties() { return surface_prop; }
227
233 std::vector<const char*>& get_surface_property_names() { return surface_prop_names; }
234};
235
242std::ifstream open_binary_file(const char* filename) {
243 return std::ifstream{filename, std::ifstream::in | std::ifstream::binary};
244}
245
253template <typename T>
254char* as_char_ptr(T* pointer) {
255 return reinterpret_cast<char*>(pointer);
256}
257
265template <typename T>
266void read_binary_value(std::ifstream& in, T* dst) {
267 in.read(as_char_ptr(dst), sizeof(T));
268}
269
278template <typename T>
279void read_binary_array(std::ifstream& in, T* dst, size_t array_length) {
280 size_t n_bytes = array_length * sizeof(T);
281 in.read(as_char_ptr(dst), n_bytes);
282}
283
289template <typename T>
290Geometry<T>::Geometry() : num_vertices(0), A_ref(1.0), c_ref(1.0), bounding_box(6, 1) {}
291
300template <typename T>
301 Geometry<T>::Geometry(std::string filename) : A_ref(1.0), c_ref(1.0) {
302
303 float p1_local[4], p2_local[4], p3_local[4], normal[4]; // Temporary arrays to hold vertex and normal data
304 char dummy[3]; // Dummy variable to skip unused bytes in the STL format
305
306 // Open the binary STL file
307 auto in = open_binary_file(filename.c_str());
308
309 // Read the STL header and the number of vertices (triangles)
310 read_binary_array<char>(in, header, 80);
311 header[80] = '\0'; // Null-terminate the header string
312 read_binary_value<unsigned int>(in, &num_vertices);
313
314 // Initialize matrices for vertices and normals based on the number of triangles
315 p1[0] = Matrix<T>(1, num_vertices);
316 p1[1] = Matrix<T>(1, num_vertices);
317 p1[2] = Matrix<T>(1, num_vertices);
318 p2[0] = Matrix<T>(1, num_vertices);
319 p2[1] = Matrix<T>(1, num_vertices);
320 p2[2] = Matrix<T>(1, num_vertices);
321 p3[0] = Matrix<T>(1, num_vertices);
322 p3[1] = Matrix<T>(1, num_vertices);
323 p3[2] = Matrix<T>(1, num_vertices);
324
325 norm_x = Matrix<T>(1, num_vertices);
326 norm_y = Matrix<T>(1, num_vertices);
327 norm_z = Matrix<T>(1, num_vertices);
328 norm_mod = Matrix<T>(1, num_vertices);
329
330 // Read the vertex and normal data for each triangle
331 for (auto i = 0; i < num_vertices; i++) {
332 read_binary_array<float>(in, normal, 3);
333 read_binary_array<float>(in, p1_local, 3);
334 read_binary_array<float>(in, p2_local, 3);
335 read_binary_array<float>(in, p3_local, 3);
336 in.read(dummy, 2); // Skip two bytes (STL format padding)
337
338 // Store the vertex and normal data in the respective matrices
339 p1[0](0, i) = static_cast<T>(p1_local[0]);
340 p1[1](0, i) = static_cast<T>(p1_local[1]);
341 p1[2](0, i) = static_cast<T>(p1_local[2]);
342 p2[0](0, i) = static_cast<T>(p2_local[0]);
343 p2[1](0, i) = static_cast<T>(p2_local[1]);
344 p2[2](0, i) = static_cast<T>(p2_local[2]);
345 p3[0](0, i) = static_cast<T>(p3_local[0]);
346 p3[1](0, i) = static_cast<T>(p3_local[1]);
347 p3[2](0, i) = static_cast<T>(p3_local[2]);
348
349 norm_x(0, i) = static_cast<T>(normal[0]);
350 norm_y(0, i) = static_cast<T>(normal[1]);
351 norm_z(0, i) = static_cast<T>(normal[2]);
352
353 // Calculate the magnitude of the normal vector and normalize it
354 norm_mod(0, i) = std::sqrt(norm_x(0, i) * norm_x(0, i) + norm_y(0, i) * norm_y(0, i) + norm_z(0, i) * norm_z(0, i));
355 norm_x(0, i) /= norm_mod(0, i);
356 norm_y(0, i) /= norm_mod(0, i);
357 norm_z(0, i) /= norm_mod(0, i);
358 }
359
360 // Compute the bounding box for the geometry
361 compute_bounding_box();
362}
363
369template <typename T>
370Geometry<T>::~Geometry() {
371 num_vertices = 0;
372}
373
380template <typename T>
381 void Geometry<T>::compute_bounding_box() {
382 this->bounding_box = Matrix<T>(6, 1);
383 T bb_x_1, bb_x_2, bb_y_1, bb_y_2, bb_z_1, bb_z_2;
384
385 // Compute the min and max values for each dimension (x, y, z)
386 bb_x_1 = std::min(p1[0].min(), std::min(p2[0].min(), p3[0].min()));
387 bb_x_2 = std::max(p1[0].max(), std::max(p2[0].max(), p3[0].max()));
388 bb_y_1 = std::min(p1[1].min(), std::min(p2[1].min(), p3[1].min()));
389 bb_y_2 = std::max(p1[1].max(), std::max(p2[1].max(), p3[1].max()));
390 bb_z_1 = std::min(p1[2].min(), std::min(p2[2].min(), p3[2].min()));
391 bb_z_2 = std::max(p1[2].max(), std::max(p2[2].max(), p3[2].max()));
392
393 // Store the bounding box values
394 bounding_box(0, 0) = bb_x_1;
395 bounding_box(1, 0) = bb_x_2;
396 bounding_box(2, 0) = bb_y_1;
397 bounding_box(3, 0) = bb_y_2;
398 bounding_box(4, 0) = bb_z_1;
399 bounding_box(5, 0) = bb_z_2;
400}
401
407template <typename T>
408 Matrix<T> Geometry<T>::get_center() {
409 Matrix<T> center(3, 1);
410
411 // Calculate the center point by averaging the min and max coordinates in each dimension
412 center(0) = (bounding_box(0, 0) + bounding_box(1, 0)) / 2.0;
413 center(1) = (bounding_box(3, 0) + bounding_box(2, 0)) / 2.0;
414 center(2) = (bounding_box(5, 0) + bounding_box(4, 0)) / 2.0;
415
416 return center;
417}
418
427template <typename T>
428 void Geometry<T>::shift(T* offset) {
429 // Shift all vertices by the specified offset in x, y, and z directions.
430 p1[0] += offset[0];
431 p1[1] += offset[1];
432 p1[2] += offset[2];
433
434 p2[0] += offset[0];
435 p2[1] += offset[1];
436 p2[2] += offset[2];
437
438 p3[0] += offset[0];
439 p3[1] += offset[1];
440 p3[2] += offset[2];
441
442 // Update bounding box by adding offset to each corresponding dimension.
443 bounding_box(0, 0) += offset[0];
444 bounding_box(1, 0) += offset[0];
445 bounding_box(2, 0) += offset[1];
446 bounding_box(3, 0) += offset[1];
447 bounding_box(4, 0) += offset[2];
448 bounding_box(5, 0) += offset[2];
449}
450
459template <typename T>
460 void Geometry<T>::shift(Matrix<T> offset) {
461 // Shift all vertices by the offset Matrix values in x, y, and z directions.
462 p1[0] += offset(0);
463 p1[1] += offset(1);
464 p1[2] += offset(2);
465
466 p2[0] += offset(0);
467 p2[1] += offset(1);
468 p2[2] += offset(2);
469
470 p3[0] += offset(0);
471 p3[1] += offset(1);
472 p3[2] += offset(2);
473
474 // Update bounding box by adding offset to each corresponding dimension.
475 bounding_box(0, 0) += offset(0);
476 bounding_box(1, 0) += offset(0);
477 bounding_box(2, 0) += offset(1);
478 bounding_box(3, 0) += offset(1);
479 bounding_box(4, 0) += offset(2);
480 bounding_box(5, 0) += offset(2);
481}
482
490template <typename T>
491 void Geometry<T>::scale(T factor) {
492 // Scale all vertices by the given factor.
493 p1[0] *= factor;
494 p1[1] *= factor;
495 p1[2] *= factor;
496
497 p2[0] *= factor;
498 p2[1] *= factor;
499 p2[2] *= factor;
500
501 p3[0] *= factor;
502 p3[1] *= factor;
503 p3[2] *= factor;
504
505 // Scale the bounding box by the given factor.
506 bounding_box *= factor;
507}
508
518template <typename T>
519 void Geometry<T>::rotate(int axis, T angle) {
520 if (axis == 0) {
521 // Rotation around the x-axis
522 Matrix<T> p1_y_new = (p1[1] * cos(angle) + p1[2] * sin(angle));
523 Matrix<T> p1_z_new = (p1[2] * cos(angle) - p1[1] * sin(angle));
524 Matrix<T> p2_y_new = (p2[1] * cos(angle) + p2[2] * sin(angle));
525 Matrix<T> p2_z_new = (p2[2] * cos(angle) - p2[1] * sin(angle));
526 Matrix<T> p3_y_new = (p3[1] * cos(angle) + p3[2] * sin(angle));
527 Matrix<T> p3_z_new = (p3[2] * cos(angle) - p3[1] * sin(angle));
528
529 Matrix<T> norm_y_new = norm_y * cos(angle) + norm_z * sin(angle);
530 Matrix<T> norm_z_new = norm_z * cos(angle) - norm_y * sin(angle);
531
532 // Update vertices and normals after the rotation.
533 p1[1] = p1_y_new; p1[2] = p1_z_new;
534 p2[1] = p2_y_new; p2[2] = p2_z_new;
535 p3[1] = p3_y_new; p3[2] = p3_z_new;
536 norm_y = norm_y_new; norm_z = norm_z_new;
537
538 compute_bounding_box(); // Recompute the bounding box after rotation.
539 } else if (axis == 1) {
540 // Rotation around the y-axis
541 Matrix<T> p1_x_new = p1[0] * cos(angle) - p1[2] * sin(angle);
542 Matrix<T> p1_z_new = p1[0] * sin(angle) + p1[2] * cos(angle);
543 Matrix<T> p2_x_new = p2[0] * cos(angle) - p2[2] * sin(angle);
544 Matrix<T> p2_z_new = p2[0] * sin(angle) + p2[2] * cos(angle);
545 Matrix<T> p3_x_new = p3[0] * cos(angle) - p3[2] * sin(angle);
546 Matrix<T> p3_z_new = p3[0] * sin(angle) + p3[2] * cos(angle);
547
548 Matrix<T> norm_x_new = (norm_x * cos(angle) - norm_z * sin(angle));
549 Matrix<T> norm_z_new = (norm_x * sin(angle) + norm_z * cos(angle));
550
551 // Update vertices and normals after the rotation.
552 p1[0] = p1_x_new; p1[2] = p1_z_new;
553 p2[0] = p2_x_new; p2[2] = p2_z_new;
554 p3[0] = p3_x_new; p3[2] = p3_z_new;
555 norm_x = norm_x_new; norm_z = norm_z_new;
556
557 compute_bounding_box(); // Recompute the bounding box after rotation.
558 } else if (axis == 2) {
559 // Rotation around the z-axis
560 Matrix<T> p1_x_new = p1[0] * cos(angle) + p1[1] * sin(angle);
561 Matrix<T> p1_y_new = p1[1] * cos(angle) - p1[0] * sin(angle);
562 Matrix<T> p2_x_new = p2[0] * cos(angle) + p2[1] * sin(angle);
563 Matrix<T> p2_y_new = p2[1] * cos(angle) - p2[0] * sin(angle);
564 Matrix<T> p3_x_new = p3[0] * cos(angle) + p3[1] * sin(angle);
565 Matrix<T> p3_y_new = p3[1] * cos(angle) - p3[0] * sin(angle);
566
567 Matrix<T> norm_x_new = norm_x * cos(angle) + norm_y * sin(angle);
568 Matrix<T> norm_y_new = norm_y * cos(angle) - norm_x * sin(angle);
569
570 // Update vertices and normals after the rotation.
571 p1[0] = p1_x_new; p1[1] = p1_y_new;
572 p2[0] = p2_x_new; p2[1] = p2_y_new;
573 p3[0] = p3_x_new; p3[1] = p3_y_new;
574 norm_x = norm_x_new; norm_y = norm_y_new;
575
576 compute_bounding_box(); // Recompute the bounding box after rotation.
577 } else {
578 std::cout << "The axis must be either 0, 1, or 2!\n";
579 }
580}
581
590template <typename T>
591T Geometry<T>::compute_area(unsigned int index) {
592 return 0.5 * fabs(p1[0](index) * (p2[1](index) * p3[2](index) - p3[1](index) * p2[2](index)) +
593 p1[1](index) * (p3[0](index) * p2[2](index) - p2[0](index) * p3[2](index)) +
594 p1[2](index) * (p2[0](index) * p3[1](index) - p3[0](index) * p2[1](index)));
595}
596
605template <typename T>
606 void Geometry<T>::add_GSI_property(const char* property_name, T value) {
607 Matrix<T> property(1, num_vertices);
608 property += value;
609 GSI_prop.push_back(property);
610 GSI_prop_names.push_back(property_name);
611}
612
621template <typename T>
622 void Geometry<T>::add_GSI_property(const char* property_name, const Matrix<T>& values) {
623 GSI_prop.push_back(values);
624 GSI_prop_names.push_back(property_name);
625}
626
635template <typename T>
636 void Geometry<T>::add_surface_property(const char* property_name, T value) {
637 Matrix<T> property(1, num_vertices);
638 property += value;
639 surface_prop.push_back(property);
640 surface_prop_names.push_back(property_name);
641}
642
651template <typename T>
652 void Geometry<T>::add_surface_property(const char* property_name, const Matrix<T>& values) {
653 surface_prop.push_back(values);
654 surface_prop_names.push_back(property_name);
655}
656
663template <typename T>
664 void Geometry<T>::print() {
665 std::cout << "STL file with " << num_vertices << " triangles:\n\n";
666
667 // Loop through each triangle and print its details
668 for (auto i = 0; i < num_vertices; i++) {
669 std::cout << "Triangle " << i << ":\n";
670 std::cout << "p1: [" << p1[0](i) << " " << p1[1](i) << " " << p1[2](i) << "]\n";
671 std::cout << "p2: [" << p2[0](i) << " " << p2[1](i) << " " << p2[2](i) << "]\n";
672 std::cout << "p3: [" << p3[0](i) << " " << p3[1](i) << " " << p3[2](i) << "]\n";
673
674 // Print the normal vector of the triangle
675 std::cout << "normal: [" << norm_x(i) << " " << norm_y(i) << " " << norm_z(i) << "]\n";
676
677 // Print the names of GSI properties
678 std::cout << "GSI properties names: ";
679 for (auto j = 0; j < GSI_prop_names.size(); j++) {
680 std::cout << GSI_prop_names[j] << " ";
681 }
682
683 // Print the values of GSI properties for the current triangle
684 std::cout << "\nGSI properties values: ";
685 for (auto j = 0; j < GSI_prop.size(); j++) {
686 std::cout << GSI_prop[j](i) << " ";
687 }
688
689 // Print the names of surface properties
690 std::cout << "\nSurface properties names: ";
691 for (auto j = 0; j < surface_prop_names.size(); j++) {
692 std::cout << surface_prop_names[j] << " ";
693 }
694
695 // Print the values of surface properties for the current triangle
696 std::cout << "\nSurface properties values: ";
697 for (auto j = 0; j < surface_prop.size(); j++) {
698 std::cout << surface_prop[j](i) << " ";
699 }
700 std::cout << "\n\n";
701 }
702}
703
711template <typename T>
712 Matrix<T> Geometry<T>::get_bounding_box() {
713 return bounding_box;
714}
Geometry
Class representing a geometric structure composed of vertices, normals, and properties.
Definition Geometry.h:11
Geometry::get_GSI_properties
std::vector< Matrix< T > > & get_GSI_properties()
Get all GSI properties.
Definition Geometry.h:212
Geometry::get_center
Matrix< T > get_center()
Get the center of the bounding box of the geometry.
Definition Geometry.h:408
Geometry::print
void print()
Print the geometry details, including vertices, normals, and properties.
Definition Geometry.h:664
Geometry::rotate
void rotate(int axis, T angle)
Rotate the geometry around a specific axis.
Definition Geometry.h:519
Geometry::shift
void shift(T *offset)
Shift the geometry by an offset.
Definition Geometry.h:428
Geometry::add_surface_property
void add_surface_property(const char *property_name, T value)
Add a scalar surface property to the geometry.
Definition Geometry.h:636
Geometry::compute_bounding_box
void compute_bounding_box()
Compute the bounding box for the geometry.
Definition Geometry.h:381
Geometry::get_surface_property_names
std::vector< const char * > & get_surface_property_names()
Get all surface property names.
Definition Geometry.h:233
Geometry::set_c_ref
void set_c_ref(T c_ref)
Set the reference speed for the geometry.
Definition Geometry.h:129
Geometry::Geometry
Geometry()
Default constructor for Geometry class.
Definition Geometry.h:290
Geometry::add_GSI_property
void add_GSI_property(const char *property_name, T value)
Add a scalar GSI property to the geometry.
Definition Geometry.h:606
Geometry::get_GSI_property
Matrix< T > & get_GSI_property(int index)
Get a GSI property by index.
Definition Geometry.h:183
Geometry::get_c_ref
T get_c_ref()
Get the reference speed.
Definition Geometry.h:148
Geometry::get_surface_properties
std::vector< Matrix< T > > & get_surface_properties()
Get all surface properties.
Definition Geometry.h:226
Geometry::get_bounding_box
Matrix< T > get_bounding_box()
Get the bounding box of the geometry.
Definition Geometry.h:712
Geometry::get_GSI_property_name
char * get_GSI_property_name(int index)
Get the name of a GSI property by index.
Definition Geometry.h:191
Geometry::get_A_ref
T get_A_ref()
Get the reference area.
Definition Geometry.h:141
Geometry::compute_area
T compute_area(unsigned int index)
Compute the area of a specific triangle in the geometry.
Definition Geometry.h:591
Geometry::set_A_ref
void set_A_ref(T A_ref)
Set the reference area for the geometry.
Definition Geometry.h:122
Geometry::scale
void scale(T factor)
Scale the geometry by a factor.
Definition Geometry.h:491
Geometry::get_GSI_property_names
std::vector< const char * > & get_GSI_property_names()
Get all GSI property names.
Definition Geometry.h:219
Geometry::~Geometry
~Geometry()
Destructor for Geometry class.
Definition Geometry.h:370
Geometry::get_num_vertices
unsigned int get_num_vertices() const
Get the number of vertices (triangles) in the geometry.
Definition Geometry.h:54
Matrix
This class implements a matrix object used for linear algebra and vectorized operations.
Definition Matrix.h:25