feat(core): New flexible states

.github/workflows/matlab.yml

@@ -20,7 +20,8 @@ jobs:
       fail-fast: false
       matrix:
         # MATLAB requires an old version of Windows
-        os: [ubuntu-latest, macos-latest]
+        # Ubuntu dropped because of libcstd++ incompatibilities
+        os: [macos-latest]
 
     steps:
       - name: Check out repository

source/Body.cpp

@@ -39,7 +39,7 @@ using namespace std;
 namespace moordyn {
 
 Body::Body(moordyn::Log* log, size_t id)
-  : io::IO(log)
+  : Instance(log)
   , bodyId(id)
 {
 #ifdef USE_VTK
@@ -273,7 +273,7 @@ Body::setDependentStates()
 		// pass above to the point and get it to calculate the forces
 		try {
 			attachedP[i]->setKinematics(rPoint, rdPoint);
-		} catch (moordyn::invalid_value_error& exception) {
+		} catch (moordyn::invalid_value_error) {
 			// Just rethrow the exception
 			throw;
 		}
@@ -313,7 +313,7 @@ Body::setDependentStates()
 		// rdRod = [baseDX, baseDY, baseDZ, endDX, endDY, endDZ]
 		try {
 			attachedR[i]->setKinematics(rRod, rdRod);
-		} catch (moordyn::invalid_value_error& exception) {
+		} catch (moordyn::invalid_value_error) {
 			// Just rethrow the exception
 			throw;
 		}
@@ -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

@@ -40,7 +40,7 @@
 #pragma once
 
 #include "Misc.hpp"
-#include "IO.hpp"
+#include "Instance.hpp"
 #include "Util/CFL.hpp"
 #include <vector>
 #include <utility>
@@ -69,11 +69,8 @@ class Rod;
  * In the configuration file the options are:
  *
  * Name/ID, X0, Y0, Z0, Xcg, Ycg, Zcg, M, V, IX, IY, IZ, CdA-x,y,z Ca-x,y,z
- *
- * moordyn::Body extends the io::IO class, allowing it to perform input/output
- * in a consistent manner.
  */
-class Body final : public io::IO, public SuperCFL
+class DECLDIR Body final : public Instance, public SuperCFL
 {
   public:
 	/** @brief Costructor
@@ -126,8 +123,6 @@ class Body final : public io::IO, 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;
 
@@ -260,15 +255,27 @@ class Body final : public io::IO, public SuperCFL
 	 */
 	inline vec6 getUnfreeVel() const { return rd_ves; }
 
-	/** @brief Initialize the FREE point state
+	/** @brief Initialize a free instance
+	 * @param r The output state variable
 	 * @return The 6-dof position (first) and the 6-dof velocity (second)
-	 * @throw moordyn::invalid_value_error If the body is not of type
-	 * moordyn::Body::FREE
+	 * @throw moordyn::invalid_value_error If the instance does not have free
+	 * states. e.g. a coupled body controlled from outside
 	 * @throws moordyn::output_file_error If an outfile has been provided, but
 	 * it cannot be written
+	 * @{
 	 */
 	std::pair<XYZQuat, vec6> initialize();
 
+	inline void initialize(InstanceStateVarView r)
+	{
+		const auto [pos, vel] = initialize();
+		r.row(0).head<7>() = pos.toVec7();
+		r.row(0).tail<6>() = vel;
+	}
+	/**
+	 * @}
+	 */
+
 	/** @brief Initialize the free body
 	 *
 	 * Those are the bodies with type moordyn::Body::FREE
@@ -378,18 +385,27 @@ class Body final : public io::IO, public SuperCFL
 	/** @brief Set the states to the body to the position r and velocity rd
 	 * @param r The position
 	 * @param rd The velocity
+	 * @{
+	 */
+	void setState(XYZQuat r, vec6 rd);
+
+	inline void setState(const InstanceStateVarView r)
+	{
+		setState(XYZQuat::fromVec7(r.row(0).head<7>()), r.row(0).tail<6>());
+	}
+	/**
+	 * @}
 	 */
-	void DECLDIR setState(XYZQuat r, vec6 rd);
 
 	/** @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
@@ -399,6 +415,18 @@ class Body final : public io::IO, 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/CMakeLists.txt

@@ -1,6 +1,7 @@
 set(MOORDYN_SRCS
     Body.cpp
     Point.cpp
+    Instance.cpp
     IO.cpp
     Line.cpp
     Log.cpp
@@ -23,6 +24,7 @@ set(MOORDYN_SRCS
 set(MOORDYN_HEADERS
     Body.hpp
     Point.hpp
+    Instance.hpp
     IO.hpp
     Line.hpp
     Log.hpp

source/IO.cpp

@@ -196,7 +196,7 @@ IO::IO(moordyn::Log* log)
   : LogUser(log)
   , _is_big_endian(false)
   , _min_major_version(2)
-  , _min_minor_version(2)
+  , _min_minor_version(4)
 {
 	_is_big_endian = is_big_endian();
 	if (_min_major_version <= MOORDYN_MAJOR_VERSION) {
@@ -206,8 +206,6 @@ IO::IO(moordyn::Log* log)
 	}
 }
 
-IO::~IO() {}
-
 void
 IO::Save(const std::string filepath)
 {

source/IO.hpp

@@ -62,7 +62,7 @@ namespace io {
  * @note Since this class already publicly inherits moordyn::LogUser, all the
  * classes inheriting this class do not require inheriting also LogUser
  */
-class IO : public LogUser
+class DECLDIR IO : public LogUser
 {
   public:
 	/** @brief Costructor
@@ -72,7 +72,7 @@ class IO : public LogUser
 
 	/** @brief Destructor
 	 */
-	~IO();
+	virtual ~IO() = default;
 
 	/** @brief Save the entity into a file
 	 *
@@ -207,6 +207,27 @@ class IO : public LogUser
 	 */
 	std::vector<uint64_t> Serialize(const std::vector<mat6>& l);
 
+	/** @brief Pack an arbitrarily large matrix
+	 * @param l The matrix
+	 * @return The packed list
+	 */
+	inline std::vector<uint64_t> Serialize(
+		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++) {
+			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
 	 * This function might act recursively
 	 * @param l The list
@@ -219,7 +240,7 @@ class IO : public LogUser
 		const uint64_t n = l.size();
 		data.push_back(Serialize(n));
 		for (auto v : l) {
-			auto subdata = Serialize(v);
+			std::vector<uint64_t> subdata = Serialize(v);
 			data.insert(data.end(), subdata.begin(), subdata.end());
 		}
 		return data;
@@ -323,6 +344,29 @@ class IO : public LogUser
 	 */
 	uint64_t* Deserialize(const uint64_t* in, std::vector<mat6>& out);
 
+	/** @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
+	 */
+	uint64_t* Deserialize(
+		const uint64_t* in,
+		Eigen::Matrix<real, Eigen::Dynamic, Eigen::Dynamic>& out)
+	{
+		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++) {
+			for (unsigned int j = 0; j < m; j++) {
+				remaining = Deserialize(remaining, out(i, j));
+			}
+		}
+		return remaining;
+	}
+
 	/** @brief Unpack a loaded list of lists
 	 *
 	 * This function might works recursively

source/Instance.cpp

@@ -0,0 +1,36 @@
+/*
+ * Copyright (c) 2019 Matt Hall <mtjhall@alumni.uvic.ca>
+ *
+ * This file is part of MoorDyn.  MoorDyn is free software: you can redistribute
+ * it and/or modify it under the terms of the GNU General Public License as
+ * published by the Free Software Foundation, either version 3 of the License,
+ * or (at your option) any later version.
+ *
+ * MoorDyn is distributed in the hope that it will be useful, but WITHOUT ANY
+ * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+
+ * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with MoorDyn.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include "Instance.hpp"
+#include <atomic>
+
+namespace moordyn {
+
+std::atomic<size_t> __instances_counter(0);
+
+Instance::Instance(moordyn::Log* log)
+  : io::IO(log)
+{
+	_id = __instances_counter++;
+}
+
+void reset_instance_ids()
+{
+	__instances_counter = 0;
+}
+
+}  // ::moordyn