I'm not hugely experienced with physics engines and am just trying to implement some water physics in Godot using Bullet physics engine, think of a boat bobbing on water. There isn't any in built explicit support for water that I am aware of, so I am trying to do it manually by applying some forces on physics ticks.

I want the boat to be able to move pretty freely in terms of yaw, but when it gets bumped I need some kind of damping on the pitch and roll to get it back to a horizontal position. As far as I can see I can read the rotation of the boat (as a 3x3 matrix or quaternion) and I can read and write the angular velocity as an xyz vector (presumably the length denotes the magnitude and the vector direction the axis of rotation in some way).

So I have 2 questions really:

1. How would be a good way of achieving this damping to pitch and roll
2. To steer the boat (yaw), if I create an angular velocity of e.g. 0, 1, 0, can I just add this to the current angular velocity?

I'm not worried about it being correct in terms of water displacement / buoyancy, just something very cheap and approximate.

2 hours ago, lawnjelly said:

(presumably the length denotes the magnitude and the vector direction the axis of rotation in some way).

The magnitude is angle per time unit (seconds, usually).

Note that with this representation you can simply add multiple angular velocities (or torques) together, so no non-commutative rotation multiplication issues here(!).

You might want to take a look on how buoyancy works: https://en.wikipedia.org/wiki/Buoyancy

I assume a typical implementation is like so: Approximate the ship with a low poly convex hull (box would do); clip this box with the water plane; calculate COM and volume of lower half and add an upwards force at the clipped COM relative to volume.

So the upforce would be this:

vec upforce = volume * someConstant * vec(0,1,0);

To apply an external force to a rigid body you calculate force and torque:

vec lever = clippedCOM - body.centerOfMassPosition; // clippedCOM is the point of application of the force

vec torque = lever.Cross(upforce); // because this point has some offset to the body COM, it causes some torque too

body.AddForce (upforce);

bodyAddTorque (torque);

2 hours ago, lawnjelly said:

To steer the boat (yaw), if I create an angular velocity of e.g. 0, 1, 0, can I just add this to the current angular velocity?

Usually you should not change velocities when using a physics engine. Imagine a box resting close to a wall and you add velocity that makes it pushing against the wall - the engine might not be able to calculate reaction forces to prevent penetration. It is better to add forces and torques instead altering velocities directly. 

But if you still want to control velocities in a simple manner without caring about physics, e.g. you want to make a PacMan game, this tool functions are very handy:

	vec ConvertLinVelToForce ( vec targetVel, vec currentVel, float timestep, float mass)
{
    vec3 force ((targetVel - currentVel) * (mass / timestep));
    return force;
}
	

So if you have a constant target velocity, this function calculates the necessary adjustment force to keep the bodies velocity constant too.

To use it you would do this:

body.AddForce( 0.3f * ConvertLinVelToForce ( pacman.targetVel, body.linearVelocity, 1/60, body.mass));

The damping factor of 0.3f is important to prevent oscillation / huge forces and the behavior depends of the value of your fixed timestep.

For the angular part the math is a bit more involved because unlike mass inertia is not uniform relative to direction:

	vec ConvertAngVelToTorque (vec targetAngVel, vec currentAngVel, matrix4x4 &matrix, float timestep, float Ixx, float Iyy, float Izz)
{
        vec torque = (targetAngVel - currentAngVel) / timestep;
        torque = matrix.Unrotate (torque); // rotate to local body space
        torque[0] *= Ixx;
        torque[1] *= Iyy; // apply nonuniform scale caused by inertia
        torque[2] *= Izz;
        torque = matrix.Rotate (torque); // rotate back to global space
        return torque;
}
	

If your engine stores inverse inertia matrix, the above could be compacted to a single matrix multiplication mostly, but this way it makes more sense i hope.

So all this is not so specific to your question, but i guess it helps a lot with 

2 hours ago, lawnjelly said:

I'm not hugely experienced with physics engines

:)

Edit: Using the 'add external force' math, you can easily model real boat motors or spaceship thrusters etc.

Blimey Joe that's a load to process I'll have to have a thorough read through in the morning.   The avenue I've been exploring in the mean time is:

• getting the the boat rotation as a quaternion
• multiplying a vector to the top of the mast in local space (0, 1, 0) by the quaternion to get the top of mast (rotated)
• using the cross product of the desired mast location (0, 1, 0) and the rotated mast to get the axis of rotation to the horizontal boat position
• using the dot product of same to get the magnitude of the rotation needed to reach this equilibrium position
• using this axis and angle to get an angular velocity to add to the current angular velocity

Still debugging to see if I can get this to work before trying something more complex.

24 minutes ago, lawnjelly said:

Blimey Joe that's a load to process I'll have to have a thorough read through in the morning.   The avenue I've been exploring in the mean time is:

• getting the the boat rotation as a quaternion
• multiplying a vector to the top of the mast in local space (0, 1, 0) by the quaternion to get the top of mast (rotated)
• using the cross product of the desired mast location (0, 1, 0) and the rotated mast to get the axis of rotation to the horizontal boat position
• using the dot product of same to get the magnitude of the rotation needed to reach this equilibrium position
• using this axis and angle to get an angular velocity to add to the current angular velocity

Still debugging to see if I can get this to work before trying something more complex.

I hear you, but real boats get upwards automatically because their com is low under the water plane, so buoyancy would result in the expected behaviour without any hacking.

But you might need the functionality to offset the COM form the geometrical center, which any physics engine provides. Or you model the boat using a compound of two bodies, a lightweight for the mast and a heavy for the bottom.

What you want to do sounds like a good opportunity to do this correctly. It's quite simple and physics is so much more fun than graphics, also gameplay should benefit a lot. (The only lame thing here seems the need to clip a polyhedron, but you could use a bunch of spheres instead.)

Okay I managed to get the 'hack' method working, it may be fine for my purposes if not I will look into Joe's method, certainly the buoyancy article was very enlightening and Archimedes method!  

For reference this is how my hack method worked:

• Intended mast position (in model space) is 0, 1, 0 (where y is up)
• Multiply ideal mast position by model rotation to get the actual current mast position
• Cross product of these two vectors to get axis of rotation to ideal position
• Angle to ideal is dot product of the 2 mast vectors, and safe acos on the result (cap dot product -1 to +1 to prevent NaN)
• Apply angular velocity towards ideal position using this axis, and the angle times by some magic value (e.g 0.1)

Although it gives easy 'righting' at the moment it pays no attention to the shape of the boat (i.e. length versus width). A real boat would presumably right on the lengthways axis quicker than on the width axis I believe, due to buoyancy? Maybe I can modify it to take account of this.

2 hours ago, lawnjelly said:

A real boat would presumably right on the lengthways axis quicker than on the width axis I believe, due to buoyancy? Maybe I can modify it to take account of this.

Yes, because the boat is more longer than wide, the clipped volume center can move forth and back more than sideways, so the torque can be larger accordingly. You could multiply the sideways dimension of your angular velocity in local ship space by 3 to fake this.

4 minutes ago, JoeJ said:

Yes, because the boat is more longer than wide, the clipped volume center can move forth and back more than sideways, so the torque can be larger accordingly. You could multiply the sideways dimension of your angular velocity in local ship space by 3 to fake this.

Sounds a good solution, thanks!  Also the magic number to be applied might need to vary with the size of the boat, I'll have to test it out.. and I could turn off the righting if the boat has capsized.

First of all define a boat hull, lets say you define it by triangles, now you cut those triangles with waterplanes to get the area that is under water, you list those triangles (newly created), then you calculate skin drag and form drag and apply that to rotation velocity, after that you need a damping function either way to cut the close to zero related problems and smoothly deaccelerate things for the simulation to run efficently.

Now pay attention

vec3 form_drag;
vec3 skin_drag;
float whole_sub_surf = 0.0;
for (int i=0; i < submerged_tri_count; i++)
	//now iterate through form_percentage_coefff
if (ACTUAL_FRAME[i].surface > 0.1) //do not coun't anything that is less than 10 cm^2
{
	whole_sub_surf = whole_sub_surf + ACTUAL_FRAME[i].surface;
	vec3 element_form_drag = vec3(0.0, 0.0, 0.0); //somehow i dont trust that compiler due to optimization will zero it
	vec3 element_skin_drag = vec3(0.0, 0.0, 0.0);
	vec3 cog2pC 	= vectorAB(pos + ROTATION_MAT * CENTER_OF_GRAVITY, ACTUAL_FRAME[i].pC);
	vec3 local_vel = vel + AngVel * cog2pC;
	vec3 vn = Normalize(local_vel);
	float vv = VectorLength(local_vel);
	float v2 = vv;
	vv = vv * vv;
float form_percentage_coeff = absnf( acos(dot(ACTUAL_FRAME[i].normal, vn)) / float(pi) );
if (dot(ACTUAL_FRAME[i].normal, vn) <= 0.0) form_percentage_coeff = 0; //suction force

element_form_drag = (-vn) * (form_percentage_coeff * (999.1026 * 0.5 * vv) * ACTUAL_FRAME[i].surface); //form drag




float skin_percentage_coeff = absnf(1.0 - form_percentage_coeff); //due to float inaccuracy
vec3 surf_n = Normal(ACTUAL_FRAME[i].V[0], ACTUAL_FRAME[i].V[1], ACTUAL_FRAME[i].V[2], true);
float surf_choordlen = getTriangleChoordLength(
-local_vel, surf_n, ACTUAL_FRAME[i].pC,
ACTUAL_FRAME[i].V[0], ACTUAL_FRAME[i].V[1], ACTUAL_FRAME[i].V[2], true, ACTUAL_FRAME[i].pC);

if (surf_choordlen > 20.0) surf_choordlen = 1.0;

//calc reynolds number
float Rn = ReynoldsNum(VectorLength(local_vel), surf_choordlen, 0.894);
bool ah = false;
if (absnf(Rn) <= 0.001 ) {Rn = 1.0; ah= true;}

float CFRn = 1.328 / (sqrt(Rn));
if (ah) CFRn = 0.0;
//float CFRn = 0.075 / (lg * lg); //1.328 / (sqrt(Rn));//

vec3 Vfi = Normalize(ProjectVectorOnPlane(ACTUAL_FRAME[i].normal, -local_vel));
element_skin_drag = (-vn) * (CFRn * ((999.1026 * vv) / 2.0 ) * ACTUAL_FRAME[i].surface * skin_percentage_coeff);



skin_drag = skin_drag + element_skin_drag;
form_drag = form_drag + element_form_drag;

cog2cob_vec = vectorAB(pos + ROTATION_MAT * CENTER_OF_GRAVITY, ACTUAL_FRAME[i].pC);

ETorque = cog2cob_vec	*  ( element_skin_drag + element_form_drag );

elements_torque = elements_torque + ETorque;

}


drag_force = form_drag + skin_drag;

result_force = gravity_force + buoyancy_force + drag_force + thrust_vec;// + rudder_force;// + drag_force + thrust_vec + rudder_force;

mat4 inv_rot = ROTATION_MAT;
inv_rot.Inverse();

vec3 local_torque = inv_rot * elements_torque;

EAngVel = EAngAcc*dt;
AngVel = AngVel + EAngVel;

vec3 vLocalAngularMoment = (AngVel*dt)*RAD_TO_DEG;



YPRangle.pitch(cos(vLocalAngularMoment.x*imopi), -sin(vLocalAngularMoment.x*imopi));
YPRangle.DoRotation();

YPRangle.yaw(cos(vLocalAngularMoment.y*imopi), sin(vLocalAngularMoment.y*imopi));
YPRangle.DoRotation();

YPRangle.roll(cos(vLocalAngularMoment.z*imopi), -sin(vLocalAngularMoment.z*imopi));
YPRangle.DoRotation();
float yawrate = 5.0*dt;
YPRangle.yaw(cos(yawrate*imopi), sin(yawrate*imopi));
YPRangle.DoRotation();

ROTATION_MAT.LoadGLMatrix(YPRangle.AIR_MATRIX);

BUO_COG = vectorAB(CENTER_OF_BUOYANCY, pos + ROTATION_MAT * CENTER_OF_GRAVITY);



vec3 accel = result_force / mass;

vel = vel + accel*dt;
pos = pos + vel*dt;

PC stands for pressure center point

epsilon = 0.001;

inline void damp( vec3 * ptr, float amount, float epsilon)
{
	(*ptr) = (*ptr) - (*ptr)/amount;
	vec3 p = vec3(absnf<float>((*ptr).x),absnf<float>((*ptr).y),absnf<float>((*ptr).z));
	if (IsNan(p.x)) p.x = 0.0;
	if (IsNan(p.y)) p.y = 0.0;
	if (IsNan(p.z)) p.z = 0.0;

	if ( (p.x > 0) && (p.x < epsilon) ) (*ptr).x = 0.0;
	if ( (p.y > 0) && (p.y < epsilon) ) (*ptr).y = 0.0;
	if ( (p.z > 0) && (p.z < epsilon) ) (*ptr).z = 0.0;
}



inline void damp( float * ptr, float amount, float epsilon)
{
	(*ptr) = (*ptr) - (*ptr)/amount;
	if (IsNan((*ptr))) (*ptr) = 0.0;
	if ( (absnf<float>(*ptr) > 0) && (absnf<float>(*ptr) < epsilon) ) (*ptr) = 0.0;
}

//---------------------------------------------------------------------------

#ifndef aeroH
#define aeroH
//---------------------------------------------------------------------------
//Mix of crappyness
//don't even think it is correct in 0.0001%
#include <math.h>
#include "DxcMath.h"

const float AIR_STD_KINEMATIC_VISCOSITY = 1.460 * pow(10.0, -5.0); // m^2/s
const float WATER_STD_KINEMATIC_VISCOSITY = 1.3083 * pow(10.0, -6.0); // m^2/s             at temp = 10*C

inline float ReynoldsNum(float V, float choord_len, float kinematic_viscosity_of_fluid)
{
return (V*choord_len) / kinematic_viscosity_of_fluid;
}

				 //0.894

inline vec3 FindTriangleAerodynamicCenter(vec3 A, vec3 B, vec3 C)
{
}

inline float getChoordLength(vec3 velocity, vec3 surf_normal, vec3 pC, vec3 * verts, int vert_count, bool project, vec3 point_on_surf)
{
//Project velocity vector onto surface
vec3 projected = Normalize(ProjectVectorOnPlane(surf_normal, velocity));

//find the biggest side and square it - this will ensure us that ray will always hit two sides
float adst = -999.0;
for (int i=0; i < vert_count; i++)
{
int next = i + 1;
if (next == vert_count) next = 0;
adst = Max(adst, n3ddistance(verts[i],verts[next]));
}
adst = adst*adst;

//compute ray
vec3 rA = pC - projected*adst;
vec3 rB = pC + projected*adst;

vec3 cp;

for (int i=0; i < vert_count; i++)
cp = cp + verts[i];

cp = cp / float(vert_count);

vec3 n  = surf_normal*1000.0;

vec3 colp;
//we compute center point to determine whenever a side face is at front or back to the center point (because i don't need to waste time and run math on paper i can just check that and use it for further processing)
int cpside;
for (int i=0; i < vert_count; i++)
{
int next = i + 1;
if (next == vert_count) next = 0;

vec3 normal = Normal(verts[i], verts[next], verts[next] + n, true);
float dst = getplaneD(normal, verts[i]);
//cpside will now determine the
cpside = classifyapointagainstaplane(cp, normal, dst); // which side faces 'the inside' of polygon  (always either front or back never on plane) - solid convex polygon
if (cpside != isBack) //when center point is not behind side plane
{
normal = -normal;
dst = getplaneD(normal, verts[i]);
}

colp = rA;
if (SegmentPlaneIntersection(normal, dst, rA, rB, colp))
{
	//if ray and face normal 'look' (face) at each other
if (dot(Normalize(vectorAB(rA,rB)), normal) < 0) rA = colp;
else rB = colp;
}

}

return n3ddistance(rA, rB);
}



inline float getTriangleChoordLength(vec3 vel, vec3 surf_n, vec3 pC, vec3 A, vec3 B, vec3 C, bool project, vec3 pop)
{
vec3 ahue[3];
ahue[0] = A;
ahue[1] = B;
ahue[2] = C;
return getChoordLength(vel, surf_n, pC, ahue, 3, project, pop);
}

inline vec3 ShearStress(vec3 Force, float Area)
{
return Force / Area;
}

//1/ (pressure*v^2*Area) * ShearStress * dot(surf_normal, tangent_surf_dir)
inline float SkinFrictionDrag( float dens, float V, float Area, float ShearStress, t3dpoint<double> surf_normal, t3dpoint<double> surf_tangent_relative_to_airflow)
{
return (dens*V*V*Area) * ShearStress * absnf(dot(surf_normal, surf_tangent_relative_to_airflow));
}

inline float BodyDrag(float V, float Area, float fluidDens, t3dpoint<double> SurfNormal, t3dpoint<double> FluidFlowDir, bool front_only)
{

			  //now im not sure if its only the percentage of the rotation or if body is rotated the area increases
			  //anyways im not increasing the area.
float ratio = -dot(FluidFlowDir, SurfNormal);
if(!front_only) ratio = absnf(ratio);


return (0.5*(V*V)*fluidDens*Area*ratio);
}

#endif

And heres a vid but keep in mind i applied rotation force and thrust force each frame so they are moving almost constantly

https://youtu.be/7YXtG73KY_k

Wow thanks wierdcat for that working solution! It does look great in your video.  I think I get the gist of how a lot of these buoyancy simulations are working, calculating the volume under water etc.

Found some good references too:

http://seawisphunter.com/blog/2015/11/24/simulating-buoyancy-part1/

http://www.randygaul.net/2014/02/14/buoyancy-rigid-bodies-and-water-surface/

https://www.gamasutra.com/view/news/237528/Water_interaction_model_for_boats_in_video_games.php

I can see how it could be quite fun calculating the buoyancy .. I might keep the hack method for this game (for now at least) as it has to run in gdscript which is quite slow, and for a large number of boats, so has to be cheap. But it looks so good I might have a go at a sailing game later with some proper calculations.