prevent body drag coeffs from being negative

source/Body.cpp

@@ -498,7 +498,7 @@ Body::doRHS()
 	cda(Eigen::seqN(0, 3)) = OrMat.transpose() * bodyCdA.head<3>();
 	cda(Eigen::seqN(3, 3)) = OrMat.transpose() * bodyCdA.tail<3>();
 	F6net +=
-	    0.5 * env->rho_w * vi.cwiseProduct(vi.cwiseAbs()).cwiseProduct(cda);
+	    0.5 * env->rho_w * vi.cwiseProduct((vi.cwiseProduct(cda)).cwiseAbs());
 
 	// Get contributions from any points attached to the body
 	for (auto attached : attachedP) {

tests/Mooring/body_tests/bodyDrag.txt

@@ -0,0 +1,30 @@
+A simple test to ensure that body drag coefficients don't go negative in certain orientations
+------------------------- LINE TYPES --------------------------------------------------
+LineType    Diam      MassDenInAir    EA        BA/-zeta    EI      tbd   Can   Cat  Cdn   Cdt  
+(-)         (m)       (kg/m)          (N)       (Pa-s/-)    (n-m^2) (-)   (-)   (-)  (-)   (-)  
+spring      0.1       10.0            500.0     -0.4        0.0     0.0   0.0   0.0  0.0   0.00
+---------------------------- BODIES -----------------------------------------------------
+ID   Attachment  X0     Y0    Z0     r0      p0     y0     Mass  CG*   I*      Volume   CdA*   Ca
+(#)     (-)      (m)    (m)   (m)   (deg)   (deg)  (deg)   (kg)  (m)  (kg-m^2)  (m^3)   (m^2)  (-)
+1       free      0     0      0     0       -90      0      1e1     0    5        0     0.1    0
+----------------------- POINTS ----------------------------------------------
+Node      Type       X         Y          Z        M        V          CdA    CA
+(-)       (-)        (m)       (m)        (m)      (kg)     (m^3)      (m^2)  (-)
+1         Body1      0.0       0          0        0        0          0      0
+2         Fixed      0.0       0         -5        0        0          0      0
+-------------------------- LINES -------------------------------------------------
+Line     LineType NodeA     NodeB  UnstrLen  NumSegs     Flags/Outputs
+(-)      (-)       (-)       (-)   (m)         (-)          (-)
+1        spring      1         2    4.5        1           -
+-------------------------- SOLVER OPTIONS---------------------------------------------------
+2        writeLog     - Write a log file
+9.81        g            - No gravity
+0.002   dtM          - time step to use in mooring integration
+3.0e6    kb           - bottom stiffness
+3.0e5    cb           - bottom damping
+70       WtrDpth      - water depth
+3.0      ICDfac       - factor by which to scale drag coefficients during dynamic relaxation IC gen
+0.015    threshIC     - threshold for IC convergence
+0.0    TmaxIC       - threshold for IC convergence
+0.01     dtIC         - Time lapse between convergence tests (s)
+--------------------------- need this line -------------------------------------------------

tests/bodies.cpp

@@ -415,6 +415,66 @@ rotatingBody(SeriesWriter* series_writer)
 
 	return true;
 }
+
+bool
+bodyDrag(SeriesWriter* series_writer)
+{
+	int err;
+	cout << endl << " => " << __PRETTY_FUNC_NAME__ << "..." << endl;
+
+	MoorDyn system = MoorDyn_Create("Mooring/body_tests/bodyDrag.txt");
+	if (!system) {
+		cerr << "Failure Creating the Mooring system" << endl;
+		return false;
+	}
+
+	unsigned int n_dof;
+	err = MoorDyn_NCoupledDOF(system, &n_dof);
+	if (err != MOORDYN_SUCCESS) {
+		cerr << "Failure getting NCoupledDOF: " << err << endl;
+		return false;
+	}
+	if (n_dof != 0) {
+		cerr << "Expected 0 DOFs but got " << n_dof << endl;
+		return false;
+	}
+
+	double f[3];
+
+	err = MoorDyn_Init(system, nullptr, nullptr);
+	if (err != MOORDYN_SUCCESS) {
+		cerr << "Failure during the mooring initialization: " << err << endl;
+		return false;
+	}
+
+	if (!write_system_vtk(system, 0, series_writer)) {
+		return false;
+	}
+
+	double t = 0, dt = 0.1;
+	double max_t = 5;
+	while (t < max_t) {
+		// do one outer time step just to make sure everything is settled
+		err = MoorDyn_Step(system, nullptr, nullptr, f, &t, &dt);
+		if (err != MOORDYN_SUCCESS) {
+			cerr << "Failure during the mooring dynamics: " << err << endl;
+			return false;
+		}
+
+		if (!write_system_vtk(system, t, series_writer)) {
+			return false;
+		}
+	}
+
+	err = MoorDyn_Close(system);
+	if (err != MOORDYN_SUCCESS) {
+		cerr << "Failure closing Moordyn: " << err << endl;
+		return false;
+	}
+
+	return true;
+}
+
 /** @brief Runs all the test
  * @return 0 if the tests have ran just fine. The index of the failing test
  * otherwise
@@ -435,6 +495,10 @@ main(int, char**)
 		return 3;
 	}
 
+	if (!bodyDrag(NULL)) {
+		return 2;
+	}
+
 	cout << "bodies.cpp passed successfully" << endl;
 	return 0;
 }