feat(core)!: Less flexible but fast extensible state variables

Author: sanguinariojoe

source/Body.cpp

@@ -472,14 +472,16 @@ Body::setState(XYZQuat pos, vec6 vel)
 	setDependentStates();
 }
 
-std::pair<XYZQuat, vec6>
-Body::getStateDeriv()
+void
+Body::getStateDeriv(InstanceStateVarView drdt)
 {
 	if ((type != FREE) && (type != CPLDPIN)) {
 		LOGERR << "getStateDeriv called for non-free body" << endl;
 		throw moordyn::invalid_value_error("Invalid body type");
 	}
 
+	XYZQuat u7;
+
 	// Get contributions from attached points (and lines attached to
 	// them) and store in FNet:
 	doRHS();
@@ -488,10 +490,10 @@ Body::getStateDeriv()
 		a6 = solveMat6(M, F6net);
 
 		// NOTE; is the above still valid even though it includes rotational DOFs?
-		dPos.pos = v6.head<3>();
+		u7.pos = v6.head<3>();
 		// this assumes that the angular velocity is about the global coordinates
 		// which is true for bodies
-		dPos.quat = 0.5 * (quaternion(0.0, v6[3], v6[4], v6[5]) * r7.quat).coeffs();
+		u7.quat = 0.5 * (quaternion(0.0, v6[3], v6[4], v6[5]) * r7.quat).coeffs();
 	} else { // Pinned body
 
 		// account for moment in response to acceleration due to inertial coupling (off-diagonal sub-matrix terms)
@@ -504,10 +506,12 @@ Body::getStateDeriv()
 		a6(Eigen::seqN(3, 3)) = M_out3.inverse() * Fnet_out3;
 
 		// dxdt = V   (velocities)
-		dPos.pos = vec::Zero();
-		dPos.quat = 0.5 * (quaternion(0.0, v6[3], v6[4], v6[5]) * r7.quat).coeffs();
+		u7.pos = vec::Zero();
+		u7.quat = 0.5 * (quaternion(0.0, v6[3], v6[4], v6[5]) * r7.quat).coeffs();
 	}
-	return std::make_pair(dPos, a6);
+
+	drdt.row(0).head<7>() = u7.toVec7();
+	drdt.row(0).tail<6>() = a6;
 };
 
 const vec6

source/Body.hpp

@@ -123,8 +123,6 @@ class DECLDIR Body final : public Instance, public SuperCFL
 	// vec6 r6;
 	/// body 6dof velocity[x/y/z]
 	vec6 v6;
-	/// body quaternion position derivative
-	XYZQuat dPos;
 	/// body 6dof acceleration[x/y/z]
 	vec6 a6;
 
@@ -379,9 +377,9 @@ class DECLDIR Body final : public Instance, public SuperCFL
 	 */
 	void setState(XYZQuat r, vec6 rd);
 
-	inline void setState(vec7 r, vec6 rd)
+	inline void setState(const InstanceStateVarView r)
 	{
-		setState(XYZQuat::fromVec7(r), rd);
+		setState(XYZQuat::fromVec7(r.row(0).head<7>()), r.row(0).tail<6>());
 	}
 	/**
 	 * @}
@@ -390,12 +388,12 @@ class DECLDIR Body final : public Instance, public SuperCFL
 	/** @brief calculate the forces and state derivatives of the body
 	 *
 	 * This function is only meant for free bodies
-	 * @return The states derivatives, i.e. the velocity (first) and the
-	 * acceleration (second)
+	 * @param drdt The states derivatives, i.e. the velocity and the
+	 * acceleration
 	 * @throw moordyn::invalid_value_error If the body is of type
 	 * moordyn::Body::FREE
 	 */
-	std::pair<XYZQuat, vec6> getStateDeriv();
+	void getStateDeriv(InstanceStateVarView drdt);
 
 	/** @brief calculates the forces on the body
 	 * @throw moordyn::invalid_value_error If the body is of type
@@ -405,6 +403,18 @@ class DECLDIR Body final : public Instance, public SuperCFL
 
 	void Output(real time);
 
+	/** @brief Get the number of state variables required by this instance
+	 * @return 1
+	 */
+	inline const size_t stateN() const { return 1; }
+
+	/** @brief Get the dimension of the state variable
+	 * @return 7 components for position quaternion and 6 components for
+	 * linear and angular velocities, i.e. 13 components
+	 * @warning This function shall be called after ::setup()
+	 */
+	inline const size_t stateDims() const { return 13; }
+
 	/** @brief Produce the packed data to be saved
 	 *
 	 * The produced data can be used afterwards to restore the saved information

source/IO.cpp

@@ -398,18 +398,6 @@ IO::Serialize(const XYZQuat& m)
 	return data;
 }
 
-std::vector<uint64_t>
-IO::Serialize(const list& m)
-{
-	std::vector<uint64_t> data;
-	const uint64_t n = m.rows();
-	data.reserve(1 + n);
-	data.push_back(Serialize(n));
-	for (uint64_t i = 0; i < n; i++)
-		data.push_back(Serialize(m(i)));
-	return data;
-}
-
 std::vector<uint64_t>
 IO::Serialize(const std::vector<real>& l)
 {
@@ -561,21 +549,6 @@ IO::Deserialize(const uint64_t* in, XYZQuat& out)
 	return remaining;
 }
 
-uint64_t*
-IO::Deserialize(const uint64_t* in, list& out)
-{
-	uint64_t n;
-	uint64_t* remaining = Deserialize(in, n);
-	if (out.rows() != n)
-		out.resize(n);
-	for (unsigned int i = 0; i < n; i++) {
-		real v;
-		remaining = Deserialize(remaining, v);
-		out(i) = v;
-	}
-	return remaining;
-}
-
 uint64_t*
 IO::Deserialize(const uint64_t* in, std::vector<real>& out)
 {

source/IO.hpp

@@ -177,12 +177,6 @@ class DECLDIR IO : public LogUser
 	 */
 	std::vector<uint64_t> Serialize(const XYZQuat& m);
 
-	/** @brief Pack a list to make it writable
-	 * @param m The list
-	 * @return The packed list
-	 */
-	std::vector<uint64_t> Serialize(const list& m);
-
 	/** @brief Pack a list of floating point numbers to make it writable
 	 * @param l The list
 	 * @return The packed list
@@ -213,23 +207,25 @@ class DECLDIR IO : public LogUser
 	 */
 	std::vector<uint64_t> Serialize(const std::vector<mat6>& l);
 
-	/** @brief Pack an arbitrarily large vector
-	 * @param l The list
+	/** @brief Pack an arbitrarily large matrix
+	 * @param l The matrix
 	 * @return The packed list
 	 */
-	template<typename T>
 	inline std::vector<uint64_t> Serialize(
-		const Eigen::Matrix<T, Eigen::Dynamic, 1>& l)
+		const Eigen::Matrix<real, Eigen::Dynamic, Eigen::Dynamic>& l)
 	{
 		std::vector<uint64_t> data;
 		const uint64_t n = l.rows();
+		const uint64_t m = l.cols();
+		data.reserve(2 + n * m);
 		data.push_back(Serialize(n));
+		data.push_back(Serialize(m));
 		for (uint64_t i = 0; i < n; i++) {
-			std::vector<uint64_t> subdata = Serialize(l(i));
-			data.insert(data.end(), subdata.begin(), subdata.end());
+			for (uint64_t j = 0; j < m; j++) {
+				data.push_back(Serialize(l(i, j)));
+			}
 		}
 		return data;
-
 	}
 
 	/** @brief Pack a list of lists to make it writable
@@ -313,13 +309,6 @@ class DECLDIR IO : public LogUser
 	 */
 	uint64_t* Deserialize(const uint64_t* in, XYZQuat& out);
 
-	/** @brief Unpack a loaded list
-	 * @param in The pointer to the next unread value
-	 * @param out The unpacked value
-	 * @return The new pointer to the remaining data to be read
-	 */
-	uint64_t* Deserialize(const uint64_t* in, list& out);
-
 	/** @brief Unpack a loaded list of floating point numbers
 	 * @param in The pointer to the next unread value
 	 * @param out The unpacked value
@@ -355,21 +344,25 @@ class DECLDIR IO : public LogUser
 	 */
 	uint64_t* Deserialize(const uint64_t* in, std::vector<mat6>& out);
 
-	/** @brief Unpack an arbitrarily large vector
+	/** @brief Unpack an arbitrarily large matrix
 	 * @param in The pointer to the next unread value
 	 * @param out The unpacked vector
 	 * @return The new pointer to the remaining data to be read
 	 */
-	template<typename Type>
-	uint64_t* Deserialize(const uint64_t* in,
-	                      Eigen::Matrix<Type, Eigen::Dynamic, 1>& out)
+	uint64_t* Deserialize(
+		const uint64_t* in,
+		Eigen::Matrix<real, Eigen::Dynamic, Eigen::Dynamic>& out)
 	{
-		uint64_t n;
-		uint64_t* remaining = Deserialize(in, n);
-		if (out.rows() != n)
-			out.resize(n);
+		uint64_t* remaining;
+		uint64_t n, m;
+		remaining = Deserialize(in, n);
+		remaining = Deserialize(remaining, m);
+		if ((out.rows() != n) || (out.cols() != m))
+			out.resize(n, m);
 		for (unsigned int i = 0; i < n; i++) {
-			remaining = Deserialize(remaining, out(i));
+			for (unsigned int j = 0; j < m; j++) {
+				remaining = Deserialize(remaining, out(i, j));
+			}
 		}
 		return remaining;
 	}

source/Instance.cpp

@@ -28,4 +28,9 @@ Instance::Instance(moordyn::Log* log)
 	_id = __instances_counter++;
 }
 
+void reset_instance_ids()
+{
+	__instances_counter = 0;
+}
+
 }  // ::moordyn

source/Instance.hpp

@@ -35,6 +35,7 @@
 
 #pragma once
 
+#include "Misc.hpp"
 #include "IO.hpp"
 #include "Log.hpp"
 #include <utility>
@@ -67,26 +68,45 @@ class DECLDIR Instance : public io::IO
 	 * The unique identifier is a growing index which identifies every single
 	 * created instance. Even if an instance is removed, the index will remain
 	 * taken.
+	 * That applies even for multiple instances generation.
 	 * @return The unique identifier
+	 * @warning It cannot be assumed that the first taken index will be 0
+	 * @warning It cannot be assumed that consequtive ids will be assigned
 	 */
-	inline const size_t Id() const { return _id; }
+	inline const size_t id() const { return _id; }
+
+	/** @brief Set the state
+	 * @param r The state variable
+	 */
+	virtual void setState(const InstanceStateVarView r) = 0;
+
+	/** @brief Calculate forces and get the derivative of the states
+	 * @param drdt Output state derivatives
+	 * @throws nan_error If nan values are detected in any node position
+	 */
+	virtual void getStateDeriv(InstanceStateVarView drdt) = 0;
 
 	/** @brief Get the number of state variables required by this instance
 	 *
 	 * This can be seen as the number of rows on the state variables holder.
 	 * @return The number of state variables
 	 */
-	virtual inline const size_t state_n() const { return 1; }
+	virtual inline const size_t stateN() const { return 1; }
 
 	/** @brief Get the dimension of the state variable.
 	 *
 	 * This can be seen as the number of columns on the state variables holder.
 	 * @return The dimension of the state variable
 	 */
-	virtual inline const size_t state_dims() const { return 6; }
+	virtual inline const size_t stateDims() const { return 6; }
   private:
 	/// Unique identifier of this instance
 	size_t _id;
 };
 
+/** @brief Reset the instances Ids, so they will be assigned again starting
+ * from 0
+ */
+void reset_instance_ids();
+
 } // ::moordyn

source/Line.cpp

@@ -679,21 +679,11 @@ Line::setState(const std::vector<vec>& pos, const std::vector<vec>& vel)
 }
 
 void
-Line::setState(const StateVarRef pos, const StateVarRef vel)
+Line::setState(const InstanceStateVarView state)
 {
-	if ((pos.rows() != N - 1) || (vel.rows() != N - 1)) {
-		LOGERR << "Invalid input size" << endl;
-		throw moordyn::invalid_value_error("Invalid input size");
-	}
 	for (unsigned int i = 1; i < N; i++) {
-		if ((pos(i - 1).rows() != 3) || (vel(i - 1).rows() != 3)) {
-			LOGERR << "Invalid point input size on line " << number
-			       << " node " << i << ". pos size is " << pos(i - 1).rows()
-			       << " and vel size is " << vel(i - 1).rows() << endl;
-			throw moordyn::invalid_value_error("Invalid input size");
-		}
-		r[i] = pos(i - 1);
-		rd[i] = vel(i - 1);
+		r[i] = state.row(i - 1).head<3>();
+		rd[i] = state.row(i - 1).segment<3>(3);
 	}
 }
 
@@ -787,7 +777,7 @@ Line::getEndSegmentMoment(EndPoints end_point, EndPoints rod_end_point) const
 }
 
 void
-Line::getStateDeriv(StateVarRef vel, StateVarRef acc)
+Line::getStateDeriv(InstanceStateVarView drdt)
 {
 	// NOTE:
 	// Jose Luis Cercos-Pita: This is by far the most consuming function of the
@@ -1200,8 +1190,8 @@ Line::getStateDeriv(StateVarRef vel, StateVarRef acc)
 		// For small systems it is usually faster to compute the inverse
 		// of the matrix. See
 		// https://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html
-		acc(i - 1) = M[i].inverse() * Fnet[i];
-		vel(i - 1) = rd[i];
+		drdt.row(i - 1).head<3>() = rd[i];
+		drdt.row(i - 1).segment<3>(3) = M[i].inverse() * Fnet[i];
 	}
 };