292 const T EPSILON = 1e-15;
293 Matrix<T> edge1(3, 1), edge2(3, 1), h(3, 1), s(3, 1), q(3, 1);
295 direction = velocity / velocity.
norm();
303 h = direction.
cross(edge2);
304 a = edge1.
tr().
dot(h)(0, 0);
307 if (a > -EPSILON && a < EPSILON) {
314 u = f * s.
tr().
dot(h)(0, 0);
316 if (u < -EPSILON || u > 1.0 + EPSILON) {
321 v = f * direction.
tr().
dot(q)(0, 0);
323 if (v < -EPSILON || u + v > 1.0 + EPSILON) {
328 T t = f * edge2.tr().dot(q)(0, 0);
333 col_position = origin + direction * t;
358 std::vector<Ray<T>> children_rays;
359 Matrix<T> p1_col(3, 1), p2_col(3, 1), p3_col(3, 1), norm_col(3, 1), pos_col(3, 1);
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();
369 T* GSI_prop_col =
new T[GSI_properties.size() + 1];
370 T* surface_prop_col =
new T[surface_properties.size() + 1];
374 unsigned int vertex_index = 0;
375 bool bounced =
false;
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)});
386 bool collided = check_vertex_collision(p1, p2, p3, norm, pos, dist);
387 if (i != spawn_vertex_index && collided && min_dist > dist) {
401 if (!bounced)
return children_rays;
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();
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);
413 T vel_in[4] = {vel_in_tan1, vel_in_tan2, vel_in_norm};
416 for (
auto i = 0; i < GSI_properties.size(); i++) {
417 GSI_prop_col[i] = GSI_properties[i](0, vertex_index);
420 for (
auto i = 0; i < surface_properties.size(); i++) {
421 surface_prop_col[i] = surface_properties[i](0, vertex_index);
425 for (
auto i = 0; i < children; i++) {
426 T vel_refl[4] = {0.0, 0.0, 0.0};
427 kernel(vel_in, vel_refl, GSI_prop_col, gas_prop_col, surface_prop_col);
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];
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);
441 *collision_index = vertex_index;
444 return children_rays;
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