Archimedean Spring

Discuss any questions about BEPUphysics or problems encountered.
Post Reply
NewMario64
Posts: 5
Joined: Wed Mar 20, 2019 3:06 pm

Archimedean Spring

Post by NewMario64 »

I am trying to add a simulation of a Archimedean spring to an older project. (The project is using BEPUphysics v1.3.0 I believe and XNA 4.0). Currently I've created a rope/chain using a series of the sphere prefabs and distant constraints (about 11 cm apart). I initialized the spheres positions to be as if they are in the wound up position of the spring. Other then position, all the spheres are Dynamic, use the default material and have the same size(about 1.5cm radius) and mass (about 200 grams I think each), except for the center most sphere is set to be a Kinematic object and never moves. finally the gravity of the space is set to (0, 980.655, 0).

What I am having trouble with is constraining the spheres so that the system is rigid enough to support it's own weight without making all the spheres imobal. Looking at the available joints, constraints etc. I believe that the correct tools are here I just don't understand how use them to make it spring. If anyone could point me at the classes I should be investigating I should be able to get it working. Also if needed I believe I can update the BEPUphysics dll.
Norbo
Site Admin
Posts: 4929
Joined: Tue Jul 04, 2006 4:45 am

Re: Archimedean Spring

Post by Norbo »

First you'll need to adjust the constraint stiffnesses. They have a default tuning which is fairly stiff for masses of about 1-10 units, so if you are using 200, you'll need to scale the stiffness and damping coefficients:

Code: Select all

            constraint.SpringSettings.Stiffness *= 10;
            constraint.SpringSettings.Damping *= 10;
Note that increasing constraint stiffness can introduce instability. The solver might not be able to converge quickly enough and the integrator might introduce excessive error.

The most powerful way of addressing this would be to change the Space.TimeStepSettings.TimeStepDuration to a smaller value and update more frequently to compensate. You can also increase the Space.Solver.IterationLimit, though that can only help so much without the shorter time steps.

(v2 tends to be more robust in these sorts of use cases, and also has a solver-substepping timestepper that can make things more efficient. Upgrading to v2 would require quite a bit of work, though- it's a completely different engine.)
NewMario64
Posts: 5
Joined: Wed Mar 20, 2019 3:06 pm

Re: Archimedean Spring

Post by NewMario64 »

Thank you for your suggestions Norbo. I was able to make my chain a bit less volatile when it falls thanks to your suggestions. specifically I've ended reducing the mass of the spheres down to 2 and I've increased the dampening on the spring settings. I have also added some static cylinders in the hope that with a little more structure I might get better results.

However I still can't seem to enforce any rigidity after the second link or so. I've been playing with a lot of different options to do so, and the most promising approach I've come across is to add a PointOnLineJoint, and a Swing Limit to my links. After their first few links the chain ends up lying on the new cylinders much like a rope might and increasing the stiffness only seems to compromise the chain's collision with the cylinders eventually allowing the chain to pass through the cylinder. (would upload pictures but I'm not sure I can yet).
I've also gotten similar results by replacing the limits and joints with WeldJoint solver groups. However this option seems to be a dead end since I can not Identify a way to strengthen the joint since it has no spring settings. Additionally, based on the classed description I imagine that the joint is not even supposed to be as mobile as it is.

(Is v2 still compatible with XNA 4.0?)
Norbo
Site Admin
Posts: 4929
Joined: Tue Jul 04, 2006 4:45 am

Re: Archimedean Spring

Post by Norbo »

However this option seems to be a dead end since I can not Identify a way to strengthen the joint since it has no spring settings. Additionally, based on the classed description I imagine that the joint is not even supposed to be as mobile as it is.
SolverGroups like the WeldJoint in v1 are just groups of other constraints, so you can adjust the spring settings of the contained BallSocketJoint and NoRotationJoint.

Could you create an isolated version in the BEPUphysicsDemos for me to look at? I could give more specific advice.
(Is v2 still compatible with XNA 4.0?)
v2 currently targets .NET Standard 2.0, which means it'll work on .NET Framework 4.7.2 (maybe 4.6.1) or .NET Core 2.0 or other standard-supporting platforms. It doesn't use any XNA types, though.
NewMario64
Posts: 5
Joined: Wed Mar 20, 2019 3:06 pm

Re: Archimedean Spring

Post by NewMario64 »

dose this work for the demo?

Main:

Code: Select all

namespace IsolatedBEPUArchimedeanSpring
{
#if WINDOWS || XBOX
    static class Program
    {
        /// <summary>
        /// The main entry point for the application.
        /// </summary>
        static void Main(string[] args)
        {
            using (ArchimedeanSpringSim game = new ArchimedeanSpringSim())
            {
                game.Run();
            }
        }
    }
#endif
}

Game:

Code: Select all

using System;
using System.Collections.Generic;
using Microsoft.Xna.Framework;
using Microsoft.Xna.Framework.Graphics;
using Microsoft.Xna.Framework.Input;
using ConversionHelper;
using BEPUphysics.Entities;
using BEPUphysics.Entities.Prefabs;
using BEPUphysics.Constraints.TwoEntity.JointLimits;
using BEPUphysics.Constraints.TwoEntity.Joints;

namespace IsolatedBEPUArchimedeanSpring
{
    public class ArchimedeanSpringSim : Microsoft.Xna.Framework.Game
    {
        private const int simPointCount = 300;

        private const float segLength = 11.0f;
        private const float cableDiameter = 3.175f;
        private const float startRadius = 67.1075f + (cableDiameter / 2.0f);

        private GraphicsDeviceManager graphics;
        private MouseCamera camera;

        private BEPUphysics.Space space;
        private List<Entity> physObjs;
        private List<Entity> colidObjs;
        private BEPUphysicsDrawer.Models.InstancedModelDrawer md;

        public ArchimedeanSpringSim()
        {
            graphics = new GraphicsDeviceManager(this);
            Content.RootDirectory = "Content";

            graphics.PreferMultiSampling = true;
            graphics.IsFullScreen = true;
            graphics.PreferredBackBufferWidth = 1920;
            graphics.PreferredBackBufferHeight = 1080;
            graphics.SynchronizeWithVerticalRetrace = true;
        }

        protected override void Initialize()
        {
            camera = new MouseCamera(graphics.GraphicsDevice.Viewport);
            List<Vector3> startPoints = CreateCoilPoints();

            physObjs = new List<Entity>();
            colidObjs = new List<Entity>();

            md = new BEPUphysicsDrawer.Models.InstancedModelDrawer(this);
            space = new BEPUphysics.Space();
            space.ForceUpdater.Gravity = new BEPUutilities.Vector3(0.0f, -980.655f, 0.0f);

            inicializePhysicsObj(startPoints);

            base.Initialize();
        }

        private List<Vector3> CreateCoilPoints()
        {
            float currentRot = 0;
            float currentWidth = 0;
            List<Vector3> startPoints = new List<Vector3>();

            startPoints.Add(new Vector3(0.0f, startRadius, 0.0f));
            double deltaRot = (MathHelper.TwoPi * segLength) / Math.Sqrt(Math.Pow(MathHelper.TwoPi * startRadius, 2.0) + Math.Pow(cableDiameter, 2.0));
            double deltaWidth = (cableDiameter * deltaRot) / MathHelper.TwoPi;
            for(int i = 1; i < simPointCount; i++)
            {
                currentWidth += (float)deltaWidth;
                currentRot += (float)deltaRot;
                startPoints.Add(new Vector3(currentWidth, (float)(Math.Cos(currentRot) * startRadius),
                                 (float)(Math.Sin(currentRot) * startRadius)));
            }
            return startPoints;
        }

        private void inicializePhysicsObj(List<Vector3> startPoints)
        {
            float cableRadius = cableDiameter / 2.0f;
            float segMass = (float)(segLength * (((Math.PI * Math.Pow(cableRadius, 2)) - (Math.PI * Math.Pow(cableRadius - 0.25, 2))) * 7.85) / 100);

            Entity prevPoint;
            BEPUutilities.Vector3 point = MathConverter.Convert(startPoints[0]);
            Entity tmpE = new Sphere(point, cableRadius);
            //tmpE.BecomeKinematic();
            prevPoint = tmpE;
            physObjs.Add(tmpE);
            space.Add(tmpE);
            md.Add(tmpE);
            for (int i = 1; i < startPoints.Count; i++)
            {
                point = MathConverter.Convert(startPoints[i]);
                tmpE = new Sphere(point, (float)cableRadius, segMass);
                physObjs.Add(tmpE);
                space.Add(tmpE);
                md.Add(tmpE);
                BEPUutilities.Vector3 axis = BEPUutilities.Vector3.Normalize(tmpE.Position - prevPoint.Position);
                Joint tmpJ = new PointOnLineJoint(prevPoint, tmpE, prevPoint.Position, axis, tmpE.Position);
                //tmpJ.SpringSettings.Stiffness *= 100;
                tmpJ.SpringSettings.Damping *= 100;
                space.Add(tmpJ);
                JointLimit tmpJL = new SwingLimit(prevPoint, tmpE, axis, axis, 0);
                //tmpJL.SpringSettings.Stiffness *= 100;
                tmpJL.SpringSettings.Damping *= 100;
                space.Add(tmpJL);
                tmpJL = new DistanceLimit(prevPoint, tmpE, prevPoint.Position, point, segLength - .01f, segLength + .01f);
                //tmpJL.SpringSettings.Stiffness *= 50.0f;
                tmpJL.SpringSettings.Damping *= 100.0f;
                space.Add(tmpJL);
                prevPoint = tmpE;
            }
            BEPUutilities.Vector3 tmpVRight = BEPUutilities.Vector3.Right;
            BEPUutilities.Vector3 tmpVUp = BEPUutilities.Vector3.Up;
            BEPUutilities.Quaternion tmpQ;
            BEPUutilities.Quaternion.GetQuaternionBetweenNormalizedVectors(ref tmpVUp, ref tmpVRight, out tmpQ);
            tmpE = new Cylinder(BEPUutilities.Vector3.Right * 12.5f, (float)25, startRadius - cableRadius);
            tmpE.Orientation = tmpQ;
            colidObjs.Add(tmpE);
            space.Add(tmpE);
            md.Add(tmpE);
            tmpE = new Cylinder(BEPUutilities.Vector3.Left * 2.5f, 5.0f, startRadius + cableDiameter);
            tmpE.Orientation = tmpQ;
            colidObjs.Add(tmpE);
            space.Add(tmpE);
            md.Add(tmpE);
            tmpE = new Cylinder(BEPUutilities.Vector3.Right * 27.5f, 5.0f, startRadius + cableDiameter);
            tmpE.Orientation = tmpQ;
            colidObjs.Add(tmpE);
            space.Add(tmpE);
            md.Add(tmpE);
        }

        protected override void Update(GameTime gameTime)
        {
            if (Keyboard.GetState().IsKeyDown(Keys.Escape))
                this.Exit();

            space.Update();
            md.Update();
            camera.Update(gameTime);

            base.Update(gameTime);
        }
        
        protected override void Draw(GameTime gameTime)
        {
            GraphicsDevice.Clear(Color.CornflowerBlue);

            md.Draw(MathConverter.Convert(camera.GetView()), MathConverter.Convert(camera.GetProjection()));

            base.Draw(gameTime);
        }
    }
}

Camera:

Code: Select all

using System;
using Microsoft.Xna.Framework;
using Microsoft.Xna.Framework.Graphics;
using Microsoft.Xna.Framework.Input;

namespace IsolatedBEPUArchimedeanSpring
{
    public class MouseCamera
    {
        private const float rotationSpeed = 0.005f;
        private const float panSpeed = 5.0f;

        private Matrix viewMatrix;
        private Matrix projectionMatrix;
        private Viewport viewPort;

        private float leftrightRot;
        private float updownRot;
        private float nearPlane;
        private float farPlane;
        private float viewAngle;
        private Vector3 cameraPosition;
        private static MouseState originalMouseState;
        private bool updateMatrix = true;

        public MouseCamera(Viewport _viewPort)
        {
            this.leftrightRot = 0;
            this.updownRot = 0;
            this.cameraPosition = new Vector3(0, 0, 200);
            viewPort = _viewPort;

            viewAngle = MathHelper.PiOver4;
            nearPlane = 1.0f;
            farPlane = 4000.0f;

            UpdateViewMatrix();

            originalMouseState = Mouse.GetState();
        }

        public void Update(GameTime gameTime)
        {
            MouseState currentMouseState = Mouse.GetState();
            if (viewPort.Bounds.Contains(currentMouseState.X, currentMouseState.Y))
            {
                if (currentMouseState.LeftButton == ButtonState.Pressed)
                {
                    float finalPanSpeed = panSpeed;

                    float xDifference = currentMouseState.X - originalMouseState.X;
                    float yDifference = currentMouseState.Y - originalMouseState.Y;
                        
                    AddToCameraPosition(new Vector3(-xDifference * finalPanSpeed, yDifference * finalPanSpeed, 0));
                    updateMatrix = true;
                }
                else if (currentMouseState.RightButton == ButtonState.Pressed)
                {
                    float xDifference = currentMouseState.X - originalMouseState.X;
                    float yDifference = currentMouseState.Y - originalMouseState.Y;
                    leftrightRot = MathHelper.WrapAngle(leftrightRot - rotationSpeed * xDifference);
                    updownRot = MathHelper.WrapAngle(updownRot - rotationSpeed * yDifference);
                    updateMatrix = true;
                }

                int scrollDiff = currentMouseState.ScrollWheelValue - originalMouseState.ScrollWheelValue;
                if (Math.Abs(scrollDiff) > 0)
                {
                    float finalPanSpeed = panSpeed;
                    AddToCameraPosition(new Vector3(0, 0, -finalPanSpeed * (float)scrollDiff));
                    updateMatrix = true;
                }
            }

            originalMouseState = currentMouseState;
            if (updateMatrix)
            {
                UpdateViewMatrix();
                updateMatrix = false;
            }
        }

        public void AddToCameraPosition(Vector3 vectorToAdd)
        {
            float moveSpeed = 0.5f;
            Matrix cameraRotation = Matrix.CreateRotationX(updownRot) * Matrix.CreateRotationY(leftrightRot);
            Vector3 rotatedVector = Vector3.Transform(vectorToAdd, cameraRotation);
            cameraPosition += moveSpeed * rotatedVector;

            UpdateViewMatrix();
        }

        public void UpdateViewMatrix()
        {
            Matrix cameraRotation = Matrix.CreateRotationX(updownRot) * Matrix.CreateRotationY(leftrightRot);

            Vector3 cameraOriginalTarget = new Vector3(0, 0, -1);
            Vector3 cameraOriginalUpVector = new Vector3(0, 1, 0);

            Vector3 cameraRotatedTarget = Vector3.Transform(cameraOriginalTarget, cameraRotation);
            Vector3 cameraFinalTarget = cameraPosition + cameraRotatedTarget;

            Vector3 cameraRotatedUpVector = Vector3.Transform(cameraOriginalUpVector, cameraRotation);
            Vector3 cameraFinalUpVector = cameraPosition + cameraRotatedUpVector;

            viewMatrix = Matrix.CreateLookAt(cameraPosition, cameraFinalTarget, cameraRotatedUpVector);
            if (nearPlane < 0)
                nearPlane = 1.0f;

            projectionMatrix = Matrix.CreatePerspectiveFieldOfView(viewAngle, viewPort.AspectRatio, nearPlane, farPlane);
        }

        public Matrix GetProjection()
        {
            return projectionMatrix;
        }

        public Matrix GetView()
        {
            return viewMatrix;
        }
    }
}

Norbo
Site Admin
Posts: 4929
Joined: Tue Jul 04, 2006 4:45 am

Re: Archimedean Spring

Post by Norbo »

This is a pathological case for iterative solvers; extremely long high stiffness dependency chains are near impossible to handle without some form of tricks in the main solvers that most game engines use. The good news is that those tricks do exist, and they can work.

The main problems are 1) getting impulses to propagate through the chain efficiently, and 2) preventing explosions caused by the integrator due to extremely high spring stiffness.

For propagation, you can create 'cheat' constraints which skip over bodies and connect more distant elements. To prevent integrator error explosions with sufficiently high stiffness, you pretty much have to use a shorter timestep.

I modified the initialize physics method to this:

Code: Select all

            List<Vector3> startPoints = CreateCoilPoints();
            Space = new Space();
            Space.ForceUpdater.Gravity = new Vector3(0.0f, -980f, 0.0f);
            //Note that we can use a lower iteration limit. Updating at an extremely high rate makes iterations less important, since the impulses don't need to travel as far in a single instant.
            Space.Solver.IterationLimit = 2;
            //The solver tries to early out sometimes, but for simulations like this, that's not a good idea- force it to finish the job.
            SolverSettings.DefaultMinimumIterationCount = Space.Solver.IterationLimit;
            //This update rate is 20x smaller than the demos default, so the demo has to take more steps to match or else everything will just run in slow motion.
            //This time step duration will allow spring frequencies of around 600hz without risking integrator explosions.
            Space.TimeStepSettings.TimeStepDuration = 1 / 1200f;

            float cableRadius = cableDiameter / 2.0f;
            float segMass = (float)(segLength * (((Math.PI * Math.Pow(cableRadius, 2)) - (Math.PI * Math.Pow(cableRadius - 0.25, 2))) * 7.85) / 100);

            Vector3 point = startPoints[0];
            Entity currentSphere = new Sphere(point, cableRadius);
            Space.Add(currentSphere);
            var previousSpheres = new QuickQueue<Entity>(new BEPUutilities.ResourceManagement.BufferPool<Entity>());
            previousSpheres.Enqueue(currentSphere);
            for (int i = 1; i < startPoints.Count; i++)
            {
                point = startPoints[i];
                currentSphere = new Sphere(point, cableRadius, segMass);
                currentSphere.LocalInertiaTensorInverse *= Matrix3x3.CreateScale(0.1f);
                Space.Add(currentSphere);
                for (int j = 0; j < previousSpheres.Count; ++j)
                {
                    var previousSphere = previousSpheres[j];
                    //Note the use of welds; they are the closest constraint to the desired behavior.
                    var weld = new WeldJoint(previousSphere, currentSphere, (previousSphere.Position + currentSphere.Position) / 2);
                    //The stiffness/damping are both scaled up to account for the gravity scale.
                    const float stiffnessScale = 10;
                    weld.BallSocketJoint.SpringSettings.Stiffness *= stiffnessScale;
                    weld.BallSocketJoint.SpringSettings.Damping *= stiffnessScale;
                    weld.NoRotationJoint.SpringSettings.Stiffness *= stiffnessScale;
                    weld.NoRotationJoint.SpringSettings.Damping *= stiffnessScale;
                    Space.Add(weld);
                }
                //Note that this allows multiple spheres to be accumulated, so a single sphere will have several constraints. These are the 'cheat' or 'skip' constraints
                //that allow impulses to propagate through the chain more easily.
                const int constraintCount = 10;
                if (previousSpheres.Count == constraintCount)
                    previousSpheres.Dequeue();
                previousSpheres.Enqueue(currentSphere);
            }
It produces a pretty nice stiff spring without any supporting cylinders.

It's also possible to handle this kind of simulation using a quality-focused solver. They tend to be much more expensive in most simulations, but for pathological cases like this, they can be faster than iterative solvers. Bepuphysics doesn't have such a solver (v1 or v2 at the moment), so these tricks are necessary. As far as game-focused engines go, I believe Newton, PhysX, presumably Havok, and possibly Bullet have some form of opt-in quality-focused solvers or subsolvers (like reduced coordinate articulations). Robotics focused packages like MuJoCo almost certainly have such solvers too. I wouldn't be much help in getting them set up, though.
Norbo
Site Admin
Posts: 4929
Joined: Tue Jul 04, 2006 4:45 am

Re: Archimedean Spring

Post by Norbo »

One thing I forgot to note in that modified version: I scaled the inertia tensor of the sphere links up by a factor of 10. It's not strictly necessary when combined with the other techniques, but it does make things easier to solve.
NewMario64
Posts: 5
Joined: Wed Mar 20, 2019 3:06 pm

Re: Archimedean Spring

Post by NewMario64 »

:D This is even better than what I was hoping for. I Already got the demo version working and I should be able to apply this in the full version fairly quickly. Thank you so much Norbo. your a bit of a life saver in this case.

I did notice the Inertia tensor but thx for mentioning it.
Last edited by NewMario64 on Tue Mar 26, 2019 8:30 pm, edited 1 time in total.
NewMario64
Posts: 5
Joined: Wed Mar 20, 2019 3:06 pm

Re: Archimedean Spring

Post by NewMario64 »

:? how do i (or do i need to) mark the thread as solved and or give props?
Norbo
Site Admin
Posts: 4929
Joined: Tue Jul 04, 2006 4:45 am

Re: Archimedean Spring

Post by Norbo »

Don't have to, glad it worked for you :)
Post Reply