uraeus/nmbd/cpp/engine/src/uraeus/numerics/geometries.cpp
// triangular_prism geometry function returning a geometry object
// --------------------------------------------------------------
triangular_prism geometry function returning a geometry object
triangular_prism geometry function returning a geometry object
geometry triangular_prism(Eigen::Vector3d& p1, Eigen::Vector3d& p2, Eigen::Vector3d& p3, double thickness)
{
    geometry geo;
    double l1, l2, l3, pr, theta, height, area, volume, density;
    Eigen::Vector3d v1 = p1 - p2 ;
    Eigen::Vector3d v2 = p1 - p3 ;
    Eigen::Vector3d v3 = p2 - p3 ;
    l1 = v1.norm();
    l2 = v2.norm();
    l3 = v3.norm();
    pr = (l1 + l2 + l3) / 2 ;
    double cos_theta = v1.transpose() * v2;
    cos_theta = cos_theta / (l1*l2) ;
    theta = std::acos(cos_theta);
    height = l2 * std::sin(theta);
    area = std::sqrt(pr*(pr - l1)*(pr-l2)*(pr-l3)) ;
    volume = area * thickness ;
    density = 7.9*1e-3;
    Eigen::Vector3d v = oriented({p1, p2, p3});
    Eigen::Matrix3d frame = triad(v, v1) ;
    double a = v2.transpose() * (v1/l1) ;
    double Ixc = (l1*std::pow(height,3)) / 36 ;
    double Iyc = ((std::pow(l1,3)*height)-(std::pow(l1,2)*height*a)+(l1*height*std::pow(a,2))) / 36 ;
    double Izc = ((std::pow(l1,3)*height)-(std::pow(l1,2)*height*a)+(l1*height*std::pow(a,2))+(l1*std::pow(height,3))) / 36 ;
    double Ix_n = (1/12) * std::pow(thickness, 3) * l1;
    double Iy_n = (1/12) * std::pow(thickness, 3) * height;
    double mass = density * volume;
    double Ix = ((mass/area) * Ixc + ((mass/thickness*l1) * Ix_n));
    double Iy = ((mass/area) * Iyc + ((mass/thickness*height) * Iy_n));
    double Iz = ((mass/area) * Izc);
    Eigen::Vector4d P = A2P(frame);
    Eigen::DiagonalMatrix<double, 3> J(Ix, Iy, Iz);
    geo.J = A(P) * J * A(P).transpose() ;
    geo.P << 1, 0, 0, 0 ;
    geo.R = centered({p1, p2, p3});
    geo.m = mass ;
    return geo;
}
geometry triangular_prism ( Eigen :: Vector3d & p1 , Eigen :: Vector3d & p2 , Eigen :: Vector3d & p3 , double thickness ) { geometry geo ; double l1 , l2 , l3 , pr , theta , height , area , volume , density ; Eigen :: Vector3d v1 = p1 - p2 ; Eigen :: Vector3d v2 = p1 - p3 ; Eigen :: Vector3d v3 = p2 - p3 ; l1 = v1 . norm ( ) ; l2 = v2 . norm ( ) ; l3 = v3 . norm ( ) ; pr = ( l1 + l2 + l3 ) / 2 ; double cos_theta = v1 . transpose ( ) * v2 ; cos_theta = cos_theta / ( l1 * l2 ) ; theta = std :: acos ( cos_theta ) ; height = l2 * std :: sin ( theta ) ; area = std :: sqrt ( pr * ( pr - l1 ) * ( pr - l2 ) * ( pr - l3 ) ) ; volume = area * thickness ; density = 7.9 * 1e-3 ; Eigen :: Vector3d v = oriented ( { p1 , p2 , p3 } ) ; Eigen :: Matrix3d frame = triad ( v , v1 ) ; double a = v2 . transpose ( ) * ( v1 / l1 ) ; double Ixc = ( l1 * std :: pow ( height , 3 ) ) / 36 ; double Iyc = ( ( std :: pow ( l1 , 3 ) * height ) - ( std :: pow ( l1 , 2 ) * height * a ) + ( l1 * height * std :: pow ( a , 2 ) ) ) / 36 ; double Izc = ( ( std :: pow ( l1 , 3 ) * height ) - ( std :: pow ( l1 , 2 ) * height * a ) + ( l1 * height * std :: pow ( a , 2 ) ) + ( l1 * std :: pow ( height , 3 ) ) ) / 36 ; double Ix_n = ( 1 / 12 ) * std :: pow ( thickness , 3 ) * l1 ; double Iy_n = ( 1 / 12 ) * std :: pow ( thickness , 3 ) * height ; double mass = density * volume ; double Ix = ( ( mass / area ) * Ixc + ( ( mass / thickness * l1 ) * Ix_n ) ) ; double Iy = ( ( mass / area ) * Iyc + ( ( mass / thickness * height ) * Iy_n ) ) ; double Iz = ( ( mass / area ) * Izc ) ; Eigen :: Vector4d P = A2P ( frame ) ; Eigen :: DiagonalMatrix < double , 3 > J ( Ix , Iy , Iz ) ; geo . J = A ( P ) * J * A ( P ) . transpose ( ) ; geo . P << 1 , 0 , 0 , 0 ; geo . R = centered ( { p1 , p2 , p3 } ) ; geo . m = mass ; return geo ; }
triangular_prism geometry function returning a geometry object
triangular_prism geometry function returning a geometry object
// The normal height of the vertix from the base. // Creating a centroidal reference frame with z-axis normal to triangle // plane and x-axis oriented with the selected base vector v1. // Calculating the principle inertia properties "moment of areas" at the // geometry centroid. // Evaluating the moments of inertia of the side-walls // Evaluating mass // Calculating the total moment of inertia from the moment of areas // Evaluate Geometry properties //  Transforming Inertia to the Global Frame