GSI_Toolbox 1.0.0
A toolbox for Gas-Surface Interaction simulations
Loading...
Searching...
No Matches
Ray.h
1
9template <typename T>
10class Ray {
11
12 private:
13 Matrix<T> origin;
14 Matrix<T> velocity;
15 T mass_flux;
16 unsigned int children;
17 int spawn_vertex_index;
18 bool failed = false;
19
20 public:
26 Ray();
27
39 Ray(Matrix<T>& origin, Matrix<T>& velocity, T m_dot, unsigned int children, int spawn_vertex_index);
40
48 Ray(const Ray<T>& other);
49
57 Ray(Ray<T>&& other);
58
64 ~Ray();
65
74 Ray& operator=(const Ray& other);
75
84 Ray& operator=(Ray&& other);
85
91 Matrix<T>& get_velocity() { return this->velocity; }
92
98 void set_velocity(Matrix<T> new_velocity) { this->velocity = new_velocity; }
99
105 Matrix<T>& get_origin() { return this->origin; }
106
112 void set_origin(Matrix<T> new_origin) { this->origin = new_origin; }
113
119 T get_mass_flux() { return mass_flux; }
120
126 void set_mass_flux(T mass_flux) { this->mass_flux = mass_flux; }
127
142 bool check_vertex_collision(Matrix<T>& p1, Matrix<T>& p2, Matrix<T>& p3, Matrix<T>& normal, Matrix<T>& col_position, T& distance);
143
151
164 std::vector<Ray<T>> propagate(Geometry<T>& geometry, const Gas<T>& gas, void (*kernel)(T*, T*, T*, T*, T*), unsigned int* collision_index);
165};
166
172template <typename T>
174 origin = Matrix<T>(3, 1, {0.0, 0.0, 0.0});
175 velocity = Matrix<T>(3, 1, {0.0, 0.0, 0.0});
176 mass_flux = 0.0;
177 children = 0;
178 spawn_vertex_index = -1;
179}
180
192template <typename T>
193Ray<T>::Ray(Matrix<T>& origin, Matrix<T>& velocity, T m_dot, unsigned int children, int spawn_vertex_index)
194 : origin(origin), velocity(velocity), mass_flux(m_dot), children(children), spawn_vertex_index(spawn_vertex_index) {}
195
203template <typename T>
204Ray<T>::Ray(const Ray<T>& other)
205 : origin(other.origin), velocity(other.velocity), children(other.children), mass_flux(other.mass_flux), spawn_vertex_index(other.spawn_vertex_index) {}
206
214template <typename T>
216 : origin(other.origin), velocity(other.velocity), children(other.children), mass_flux(other.mass_flux), spawn_vertex_index(other.spawn_vertex_index) {
217 other.children = 0;
218 other.mass_flux = 0.0;
219 other.spawn_vertex_index = 0;
220}
221
227template <typename T>
229 children = 0;
230 mass_flux = 0.0;
231}
232
241template <typename T>
243 origin = other.origin;
244 velocity = other.velocity;
245 mass_flux = other.mass_flux;
246 children = other.children;
247 spawn_vertex_index = other.spawn_vertex_index;
248 return *this;
249}
250
259template <typename T>
261 origin = other.origin;
262 velocity = other.velocity;
263 mass_flux = other.mass_flux;
264 children = other.children;
265 spawn_vertex_index = other.spawn_vertex_index;
266
267 // Reset the properties of the moved ray
268 other.children = 0;
269 other.mass_flux = 0.0;
270 other.spawn_vertex_index = 0;
271 return *this;
272}
273
289template <typename T>
290bool Ray<T>::check_vertex_collision(Matrix<T>& p1, Matrix<T>& p2, Matrix<T>& p3, Matrix<T>& normal, Matrix<T>& col_position, T& distance) {
291
292 const T EPSILON = 1e-15; // Small tolerance for floating-point comparisons
293 Matrix<T> edge1(3, 1), edge2(3, 1), h(3, 1), s(3, 1), q(3, 1);
294 Matrix<T> direction(3, 1);
295 direction = velocity / velocity.norm();
296 T a, f, u, v;
297
298 // Compute the two edge vectors of the triangle
299 edge1 = p2 - p1;
300 edge2 = p3 - p1;
301
302 // Calculate the determinant (cross product of ray direction and edge2)
303 h = direction.cross(edge2);
304 a = edge1.tr().dot(h)(0, 0);
305
306 // Check if the ray is parallel to the triangle
307 if (a > -EPSILON && a < EPSILON) {
308 return false; // Parallel ray
309 }
310
311 // Calculate the ray-triangle intersection using barycentric coordinates
312 f = 1.0 / a;
313 s = (origin - p1);
314 u = f * s.tr().dot(h)(0, 0);
315
316 if (u < -EPSILON || u > 1.0 + EPSILON) {
317 return false; // No intersection
318 }
319
320 q = s.cross(edge1);
321 v = f * direction.tr().dot(q)(0, 0);
322
323 if (v < -EPSILON || u + v > 1.0 + EPSILON) {
324 return false; // No intersection
325 }
326
327 // Calculate the intersection point
328 T t = f * edge2.tr().dot(q)(0, 0);
329
330 if (t > EPSILON) {
331 // Ray intersects the triangle
332 distance = fabs(t);
333 col_position = origin + direction * t;
334 return true;
335 } else {
336 return false; // Intersection occurs, but not in the direction of the ray
337 }
338}
339
355template <typename T>
356std::vector<Ray<T>> Ray<T>::propagate(Geometry<T>& geometry, const Gas<T>& gas, void (*kernel)(T*, T*, T*, T*, T*), unsigned int* collision_index) {
357
358 std::vector<Ray<T>> children_rays; // Vector to store the child rays after propagation
359 Matrix<T> p1_col(3, 1), p2_col(3, 1), p3_col(3, 1), norm_col(3, 1), pos_col(3, 1); // Collision-related matrices
360 Matrix<T> p1_x = geometry.get_p1_x(), p1_y = geometry.get_p1_y(), p1_z = geometry.get_p1_z();
361 Matrix<T> p2_x = geometry.get_p2_x(), p2_y = geometry.get_p2_y(), p2_z = geometry.get_p2_z();
362 Matrix<T> p3_x = geometry.get_p3_x(), p3_y = geometry.get_p3_y(), p3_z = geometry.get_p3_z();
363 Matrix<T> norm_x = geometry.get_norm_x(), norm_y = geometry.get_norm_y(), norm_z = geometry.get_norm_z();
364 std::vector<Matrix<T>> GSI_properties(geometry.get_GSI_properties()); // GSI (General Surface Interaction) properties
365 std::vector<const char*> GSI_property_names(geometry.get_GSI_property_names()); // Names of GSI properties
366 std::vector<Matrix<T>> surface_properties(geometry.get_surface_properties()); // Surface properties at each vertex
367 std::vector<const char*> surface_property_names(geometry.get_surface_property_names()); // Names of surface properties
368
369 T* GSI_prop_col = new T[GSI_properties.size() + 1]; // GSI properties at the collision point
370 T* surface_prop_col = new T[surface_properties.size() + 1]; // Surface properties at the collision point
371 T gas_prop_col[5] = {gas.get_temperature(), gas.get_density(), gas.get_molar_mass(), gas.get_speed()}; // Gas properties
372
373 T min_dist = 1e15; // Initialize minimum distance with a large value
374 unsigned int vertex_index = 0; // Vertex index where the collision occurred
375 bool bounced = false; // Flag indicating whether the ray has collided with a surface
376
377 // Loop through each vertex in the geometry to check for collisions
378 for (auto i = 0; i < geometry.get_num_vertices(); i++) {
379 Matrix<T> p1(3, 1, {p1_x(i), p1_y(i), p1_z(i)});
380 Matrix<T> p2(3, 1, {p2_x(i), p2_y(i), p2_z(i)});
381 Matrix<T> p3(3, 1, {p3_x(i), p3_y(i), p3_z(i)});
382 Matrix<T> norm(3, 1, {norm_x(i), norm_y(i), norm_z(i)});
383 Matrix<T> pos(3, 1);
384
385 T dist = 0.0;
386 bool collided = check_vertex_collision(p1, p2, p3, norm, pos, dist); // Check if the ray intersects the triangle
387 if (i != spawn_vertex_index && collided && min_dist > dist) {
388 // Update the collision data if this collision is closer than previous ones
389 p1_col = p1;
390 p2_col = p2;
391 p3_col = p3;
392 norm_col = norm;
393 pos_col = pos;
394 vertex_index = i;
395 min_dist = dist;
396 bounced = true; // Set the flag indicating a collision has occurred
397 }
398 }
399
400 // If no collision occurred, return an empty vector (no child rays generated)
401 if (!bounced) return children_rays;
402
403 // Compute the reflection directions after collision
404 Matrix<T> norm_direction = norm_col;
405 Matrix<T> tan1_direction = (velocity - norm_direction * velocity.tr().dot(norm_direction)(0, 0)) / (velocity - norm_direction * velocity.tr().dot(norm_direction)(0, 0)).norm();
406 Matrix<T> tan2_direction = norm_col.cross(tan1_direction); // Tangential direction orthogonal to the normal
407
408 // Compute the velocity components along the normal and tangential directions
409 T vel_in_norm = (velocity.tr().dot(norm_direction))(0, 0);
410 T vel_in_tan1 = (velocity.tr().dot(tan1_direction))(0, 0);
411 T vel_in_tan2 = (velocity.tr().dot(tan2_direction))(0, 0);
412
413 T vel_in[4] = {vel_in_tan1, vel_in_tan2, vel_in_norm}; // Store input velocities for kernel function
414
415 // Extract the GSI and surface properties at the point of collision
416 for (auto i = 0; i < GSI_properties.size(); i++) {
417 GSI_prop_col[i] = GSI_properties[i](0, vertex_index);
418 }
419
420 for (auto i = 0; i < surface_properties.size(); i++) {
421 surface_prop_col[i] = surface_properties[i](0, vertex_index);
422 }
423
424 // Generate child rays after collision
425 for (auto i = 0; i < children; i++) {
426 T vel_refl[4] = {0.0, 0.0, 0.0}; // Reflected velocity components
427 kernel(vel_in, vel_refl, GSI_prop_col, gas_prop_col, surface_prop_col); // Apply the kernel function to compute reflection
428
429 // Compute the reflected velocity vector for the child ray
430 Matrix<T> child_velocity(3, 1);
431 child_velocity = norm_direction * (-vel_refl[2]) * sgn<T>(vel_in[2]) +
432 tan1_direction * vel_refl[0] +
433 tan2_direction * vel_refl[1];
434
435 // Create the child ray and add it to the list of child rays
436 Ray<T> child_ray(pos_col, child_velocity, mass_flux / children * velocity.norm() / child_velocity.norm(), 1, vertex_index);
437 children_rays.push_back(child_ray);
438 }
439
440 // Store the index of the vertex where the collision occurred
441 *collision_index = vertex_index;
442
443 // Return the generated child rays
444 return children_rays;
445}
Class representing a Gas with properties such as density, temperature, molar mass,...
Definition Gas.h:10
T get_molar_mass() const
Returns the molar mass of the gas.
Definition Gas.h:141
T get_density() const
Returns the density of the gas.
Definition Gas.h:134
T get_temperature() const
Returns the temperature of the gas.
Definition Gas.h:127
T get_speed() const
Returns the speed of the gas.
Definition Gas.h:148
Class representing a geometric structure composed of vertices, normals, and properties.
Definition Geometry.h:11
std::vector< Matrix< T > > & get_GSI_properties()
Get all GSI properties.
Definition Geometry.h:212
std::vector< const char * > & get_surface_property_names()
Get all surface property names.
Definition Geometry.h:233
std::vector< Matrix< T > > & get_surface_properties()
Get all surface properties.
Definition Geometry.h:226
std::vector< const char * > & get_GSI_property_names()
Get all GSI property names.
Definition Geometry.h:219
unsigned int get_num_vertices() const
Get the number of vertices (triangles) in the geometry.
Definition Geometry.h:54
This class implements a matrix object used for linear algebra and vectorized operations.
Definition Matrix.h:25
Matrix cross(const Matrix &other)
Calculates the cross product of two matrices.
Definition Matrix.h:624
Matrix dot(const Matrix &other)
Calculates the dot product of two matrices.
Definition Matrix.h:604
T norm()
Calculates the Frobenius norm of the matrix.
Definition Matrix.h:674
Matrix tr()
Transposes the matrix.
Definition Matrix.h:705
Class representing a Ray for simulating particle movement in a geometric domain.
Definition Ray.h:10
void set_mass_flux(T mass_flux)
Sets the mass flux of the ray.
Definition Ray.h:126
void set_origin(Matrix< T > new_origin)
Sets the origin of the ray.
Definition Ray.h:112
Matrix< T > & get_velocity()
Gets the velocity of the ray.
Definition Ray.h:91
T get_mass_flux()
Gets the mass flux of the ray.
Definition Ray.h:119
bool check_vertex_collision(Matrix< T > &p1, Matrix< T > &p2, Matrix< T > &p3, Matrix< T > &normal, Matrix< T > &col_position, T &distance)
Checks if the ray collides with a triangle defined by three vertices.
Definition Ray.h:290
void set_velocity(Matrix< T > new_velocity)
Sets the velocity of the ray.
Definition Ray.h:98
Ray()
Default constructor for the Ray class.
Definition Ray.h:173
Ray & operator=(const Ray &other)
Copy assignment operator for the Ray class.
Definition Ray.h:242
std::vector< Ray< T > > propagate(Geometry< T > &geometry, const Gas< T > &gas, void(*kernel)(T *, T *, T *, T *, T *), unsigned int *collision_index)
Propagates the ray through the geometry and handles collision and reflection.
Definition Ray.h:356
bool check_bb_collision(Matrix< T > &bb)
Checks if the ray collides with a bounding box.
~Ray()
Destructor for the Ray class.
Definition Ray.h:228
Matrix< T > & get_origin()
Gets the origin of the ray.
Definition Ray.h:105