| Re: [eigen] Optimization advice for a specific expression |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]
I am reattaching the code, as it did not get sent correctly:
#include <iostream>
#include <ctime>
#include <chrono>
#include <Eigen/Core>
#include <Eigen/Geometry>
// @InProceedings{Euler2006,
// Title = {Fast Computation Of Inertia Through
// Affinely Extended Euler And Tensor},
// Author = {DiCarlo, Antonio and Paoluzzi, Alberto},
// Year = {2006}}
template<class REAL, class INT>
class InertiaTensor
{
public:
// ----------------------------------------------------------------------------------------- Constructor
/// Constructor without parameters.
InertiaTensor()
: dens_(0.0)
{
Eplus << 2 , 1 , 1 , 5,
1 , 2 , 1 , 5,
1 , 1 , 2 , 5,
5 , 5 , 5 , 20; // 120 * E+
}
// ------------------------------------------------------------------------------------------- Accessors
/// Set the density of the volume.
/// \param dens Density of the volume
inline void setDensity(const REAL dens){ dens_ = dens; };
/// Returns the mass of the volume
/// \return The volume mass, a scalar
inline REAL getMass() { return mass_; };
/// Returns the position of the centre of mass of the volume (defined in the same coordinates as the mesh)
/// \return A 3x1 array
inline REAL* getrG() { return rG_; };
/// Returns the inertia tensor of the volume about its centre of mass (defined in the same coordinates as the mesh)
/// \return A 3x3 array
inline REAL* getInertiaTensor() { return J_; };
// ------------------------------------------------------------------------------------------- Begin Vol
/// Starts the evaluation of the volumetric properties of a new volume.
/// Clears the integrals, mass, centre of mass position, and inertia tensor.
void beginVolume()
{
InertiaTensor4x4.setZero();
}
// --------------------------------------------------------------------------------------------- End Vol
/// Finalizes the volume integration and evaluates the physical properties.
void endVolume()
{
computeMassProperties();
}
// --------------------------------------------------------------------------------------------- AddFace
/// Adds a new face to the integration process.
/// \param vertex_data A 3xnum_vertices array with the x, y, and z mesh coordinates of the face vertices.
/// \param normal_data A 3x1 array with the x, y, and z mesh coordinates of the face normal. May be set to NULL if the face is a triangle.
/// \param num_vertices Number of vertices in the face.
void addFace(const REAL* vertex_data, const REAL* normal_data, const INT num_vertices)
{
// Check that we have at least 2 vertices per face
assert(num_vertices == 3);
Eigen::Map<const Eigen::Matrix<REAL, 3, 3>> v(vertex_data);
Eigen::Matrix<REAL, 4, 4> G;
// G = [v0 v1 v2 0]
// [ 0 0 0 1]
G.col(0).head(3) = v.col(0);
G.col(1).head(3) = v.col(1);
G.col(2).head(3) = v.col(2);
G.col(3).setZero();
G.row(3).setZero();
G(3, 3) = 1.0;
Eigen::Matrix<REAL, 4, 4> temp;
Eigen::Vector4d w = G.template leftCols<3>().rowwise().sum();
temp << G.template block<3, 3>(0, 0) * G.template block<3, 3>(0, 0).transpose() + w.template head<3>() * w.transpose().template head<3>(), 5 * w.template head<3>(),
5*w.transpose().template head<3>(), 20;
// I4x4 = det(G) G E+ Gᵀ
InertiaTensor4x4.noalias() += temp * (G.template block<3,3>(0,0).determinant() / 120.0);
// Initial implementation
// InertiaTensor4x4.noalias() += (G.template block<3,3>(0,0).determinant() / 120.0) * G * Eplus * G.transpose();
}
private:
REAL dens_; // Density
REAL mass_; // Mass of the volume
REAL rG_[3]; // Position of the centre of mass
REAL J_[9]; // Inertia tensor
Eigen::Matrix<REAL, 4, 4> Eplus, InertiaTensor4x4;
// ----------------------------------------------------------------------------- Compute mass properties
void computeMassProperties()
{
// Evaluate mass
mass_ = InertiaTensor4x4(3, 3) * dens_;
}
};
template<class F>
void massProps(const char *object, const char* method, unsigned int niters, double dens, double *mass1, double *J1, double *rg1, int nf, const F &feval) {
std::cout << "--------------------------------------------------" << std::endl;
std::cout << object << ": " << nf << " faces." << std::endl;
std::cout << "--------------------------------------------------" << std::endl;
std::chrono::high_resolution_clock::time_point t1 = std::chrono::high_resolution_clock::now();
for (unsigned int iters = 0; iters < niters; iters++)
{
feval();
}
std::chrono::high_resolution_clock::time_point t2 = std::chrono::high_resolution_clock::now();
std::chrono::duration<double> time_span1 = std::chrono::duration_cast<std::chrono::duration<double>>(t2 - t1);
std::cout << object << " - " << method << std::endl;
std::cout << "Elapsed time: " << time_span1.count() << " seconds." << std::endl;
std::cout << "Total mass: " << *mass1 << std::endl;
std::cout << std::endl;
}
template <class DT>
int computeMassPropertiesIT(DT *v, int *i, int nf, const double dens,
double *mass, double *r, double *J)
{
InertiaTensor<double, int> vp;
vp.setDensity(dens);
double faceverts[9];
Eigen::Map<Eigen::Matrix<DT, 3 , Eigen::Dynamic>> verts(v, 3, *std::max_element(i, i + nf) + 1);
vp.beginVolume();
for (int face = 0; face < nf; face++)
{
Eigen::Map<Eigen::Matrix3d> points(faceverts);
points.col(0) = verts.col(i[face * 3 + 0]).template cast<double>();
points.col(1) = verts.col(i[face * 3 + 1]).template cast<double>();
points.col(2) = verts.col(i[face * 3 + 2]).template cast<double>();
vp.addFace(faceverts, nullptr, 3);
}
vp.endVolume();
*mass = vp.getMass();
Eigen::Map<Eigen::Vector3d>(r, 3, 1) = Eigen::Map<Eigen::Vector3d>(vp..getrG(), 3, 1);
Eigen::Map<Eigen::Matrix3d>(J, 3, 3) = Eigen::Map<Eigen::Matrix3d>(vp..getInertiaTensor(), 3, 3);
}
int main(int argc, char *argv[])
{
// Definition of a tetrahedron ____________________________________________________________
// Number of faces
int nf = 4;
// Array of normals
double n[12] = {
1.0 / sqrt(3.0), 1.0 / sqrt(3.0), 1.0 / sqrt(3.0),
0.0, -1.0, 0.0,
0.0, 0.0, -1.0,
-1.0, 0.0, 0.0,
};
// Array of vertices (x, y, z coordinates of each vertex)
double v[12] = {
0.0, 0.0, 0.0,
1.0, 0.0, 0.0,
0.0, 1.0, 0.0,
0.0, 0.0, 1.0
};
// Array of vertex indices (vertices in each face)
int i[12] = {
1, 2, 3,
0, 1, 3,
0, 2, 1,
0, 3, 2
};
// Evaluation of mass properties __________________________________________________________
double rg1[3], J1[9];
double mass1;
const double dens = 7780.0;
const unsigned int niters_pyramid = 1000000;
std::cout << "Using Eigen version: " << EIGEN_WORLD_VERSION << "." << EIGEN_MAJOR_VERSION << "." << EIGEN_MINOR_VERSION << "\n";
massProps("PYRAMID", "DICARLO", niters_pyramid, dens, &mass1, J1, rg1, nf,
[&](){ computeMassPropertiesIT(v, i, nf, dens, &mass1, rg1, J1);});
return 0;
}