// // This unit is part of the GLScene Engine, http://glscene.org // unit GLVerletTypes; (* Base Verlet modelling/simulation classes. This unit is generic, GLScene-specific sub-classes are in GLVerletClasses. Note that currently, the SatisfyConstraintForEdge methods push the nodes in the edge uniformly - it should push the closer node more for correct physics. It's a matter of leverage. *) interface {$I GLScene.inc} uses System.Classes, System.SysUtils, System.Types, GLObjects, GLScene, GLCoordinates, GLVectorGeometry, GLVectorLists, GLSpacePartition, GLGeometryBB, GLVectorTypes; const G_DRAG = 0.0001; cDEFAULT_CONSTRAINT_FRICTION = 0.6; type TGLVerletEdgeList = class; TGLVerletWorld = class; TVerletProgressTimes = packed record deltaTime, newTime: Double; sqrDeltaTime, invSqrDeltaTime: Single; end; // Basic verlet node TVerletNode = class(TSpacePartitionLeaf) private FForce: TAffineVector; FOwner: TGLVerletWorld; FWeight, FInvWeight: Single; FRadius: Single; FNailedDown: Boolean; FFriction: Single; FChangedOnStep: Integer; function GetSpeed: TAffineVector; protected FLocation, FOldLocation: TAffineVector; procedure SetLocation(const Value: TAffineVector); virtual; procedure SetWeight(const Value: Single); procedure AfterProgress; virtual; public constructor CreateOwned(const aOwner: TGLVerletWorld); virtual; destructor Destroy; override; // Applies friction procedure ApplyFriction(const friction, penetrationDepth: Single; const surfaceNormal: TAffineVector); // Simple and less accurate method for friction procedure OldApplyFriction(const friction, penetrationDepth: Single); // Perform Verlet integration procedure Verlet(const vpt: TVerletProgressTimes); virtual; (* Initlializes the node. For the base class, it just makes sure that FOldPosition = FPosition, so that speed is zero *) procedure Initialize; dynamic; // Calculates the distance to another node function DistanceToNode(const node: TVerletNode): Single; // Calculates the movement of the node function GetMovement: TAffineVector; (* The TVerletNode inherits from TSpacePartitionLeaf, and it needs to know how to publish itself. The owner ( a TGLVerletWorld ) has a spatial partitioning object *) procedure UpdateCachedAABBAndBSphere; override; { The VerletWorld that owns this verlet } property Owner: TGLVerletWorld read FOwner; { The location of the verlet } property Location: TAffineVector read FLocation write SetLocation; { The old location of the verlet. This is used for verlet integration } property OldLocation: TAffineVector read FOldLocation write FOldLocation; { The radius of the verlet node - this has been more or less deprecated } property Radius: Single read FRadius write FRadius; { A sum of all forces that has been applied to this verlet node during a step } property Force: TAffineVector read FForce write FForce; { If the node is nailed down, it can't be moved by either force, constraint or verlet integration - but you can still move it by hand } property NailedDown: Boolean read FNailedDown write FNailedDown; { The weight of a node determines how much it's affected by a force } property Weight: Single read FWeight write SetWeight; { InvWeight is 1/Weight, and is kept up to date automatically } property InvWeight: Single read FInvWeight; { Returns the speed of the verlet node. Speed = Movement / deltatime } property Speed: TAffineVector read GetSpeed; { Each node has a friction that effects how it reacts during contacts. } property friction: Single read FFriction write FFriction; { What phyisics step was this node last changed? Used to keep track of when the spatial partitioning needs to be updated } property ChangedOnStep: Integer read FChangedOnStep; end; TGLVerletNodeClass = class of TVerletNode; TVerletNodeList = class(TList) private function GetItems(i: Integer): TVerletNode; procedure SetItems(i: Integer; const Value: TVerletNode); public property Items[i: Integer]: TVerletNode read GetItems write SetItems; default; end; TVerletConstraint = class(TObject) private FOwner: TGLVerletWorld; FEnabled: Boolean; FTag: Integer; public constructor Create(const aOwner: TGLVerletWorld); virtual; destructor Destroy; override; { Updates the position of one or several nodes to make sure that they don't violate the constraint } procedure SatisfyConstraint(const iteration, maxIterations: Integer); virtual; abstract; { Notifies removal of a node } procedure RemoveNode(const aNode: TVerletNode); virtual; abstract; { Method that's fired before the physics iterations are performed } procedure BeforeIterations; virtual; { Onwer of the constraint } property Owner: TGLVerletWorld read FOwner; { Determines if the constraint should be enforced or not } property Enabled: Boolean read FEnabled write FEnabled; { Tag field reserved for the user. } property Tag: Integer read FTag write FTag; end; TGLVerletDualConstraint = class(TVerletConstraint) private FNodeA, FNodeB: TVerletNode; public procedure RemoveNode(const aNode: TVerletNode); override; { Reference to NodeA. } property NodeA: TVerletNode read FNodeA write FNodeA; { Reference to NodeB. } property NodeB: TVerletNode read FNodeB write FNodeB; end; TVerletGroupConstraint = class(TVerletConstraint) private FNodes: TVerletNodeList; public constructor Create(const aOwner: TGLVerletWorld); override; destructor Destroy; override; procedure RemoveNode(const aNode: TVerletNode); override; { The list of nodes that this constraint will effect } property Nodes: TVerletNodeList read FNodes; end; // Verlet edges simulate rigid collission edges TGLVerletEdge = class(TSpacePartitionLeaf) private FNodeA: TVerletNode; FNodeB: TVerletNode; public { The TGLVerletEdge inherits from TSpacePartitionLeaf, and it needs to know how to publish itself. The owner ( a TGLVerletWorld ) has a spatial partitioning object } procedure UpdateCachedAABBAndBSphere; override; constructor CreateEdgeOwned(const aNodeA, aNodeB: TVerletNode); { One of the nodes in the edge } property NodeA: TVerletNode read FNodeA write FNodeA; { One of the nodes in the edge } property NodeB: TVerletNode read FNodeB write FNodeB; end; TGLVerletEdgeList = class(TList) private function GetItems(i: Integer): TGLVerletEdge; procedure SetItems(i: Integer; const Value: TGLVerletEdge); public property Items[i: Integer]: TGLVerletEdge read GetItems write SetItems; default; end; TGLVerletGlobalConstraint = class(TVerletConstraint) private FKickbackForce: TAffineVector; FKickbackTorque: TAffineVector; FLocation: TAffineVector; procedure SetLocation(const Value: TAffineVector); virtual; public constructor Create(const aOwner: TGLVerletWorld); override; destructor Destroy; override; procedure RemoveNode(const aNode: TVerletNode); override; procedure BeforeIterations; override; procedure SatisfyConstraint(const iteration, maxIterations: Integer); override; procedure SatisfyConstraintForNode(const aNode: TVerletNode; const iteration, maxIterations: Integer); virtual; abstract; procedure SatisfyConstraintForEdge(const aEdge: TGLVerletEdge; const iteration, maxIterations: Integer); virtual; property Location: TAffineVector read FLocation write SetLocation; { The force that this collider has experienced while correcting the verlet possitions. This force can be applied to ODE bodies, for instance } property KickbackForce: TAffineVector read FKickbackForce write FKickbackForce; { The torque that this collider has experienced while correcting the verlet possitions, in reference to the center of the collider. The torque force can be applied to ODE bodies, but it must first be translated. A torque can be trasnalted by EM(b) = EM(a) + EF x VectorSubtract(b, a). Simply adding the torque to the body will NOT work correctly. See TranslateKickbackTorque } property KickbackTorque: TAffineVector read FKickbackTorque write FKickbackTorque; procedure AddKickbackForceAt(const Pos: TAffineVector; const Force: TAffineVector); function TranslateKickbackTorque(const TorqueCenter: TAffineVector): TAffineVector; end; TGLVerletGlobalFrictionConstraint = class(TGLVerletGlobalConstraint) private FFrictionRatio: Single; public constructor Create(const aOwner: TGLVerletWorld); override; property FrictionRatio: Single read FFrictionRatio write FFrictionRatio; end; TGLVerletGlobalFrictionConstraintSP = class(TGLVerletGlobalFrictionConstraint) public procedure SatisfyConstraint(const iteration, maxIterations: Integer); override; procedure PerformSpaceQuery; virtual; abstract; end; TGLVerletGlobalFrictionConstraintSphere = class(TGLVerletGlobalFrictionConstraintSP) private FCachedBSphere: TBSphere; procedure SetLocation(const Value: TAffineVector); override; public procedure UpdateCachedBSphere; procedure PerformSpaceQuery; override; function GetBSphere: TBSphere; virtual; abstract; property CachedBSphere: TBSphere read FCachedBSphere; end; TGLVerletGlobalFrictionConstraintBox = class(TGLVerletGlobalFrictionConstraintSP) private FCachedAABB: TAABB; procedure SetLocation(const Value: TAffineVector); override; public procedure UpdateCachedAABB; procedure PerformSpaceQuery; override; function GetAABB: TAABB; virtual; abstract; property CachedAABB: TAABB read FCachedAABB; end; TVerletConstraintList = class(TList) private function GetItems(i: Integer): TVerletConstraint; procedure SetItems(i: Integer; const Value: TVerletConstraint); public property Items[i: Integer]: TVerletConstraint read GetItems write SetItems; default; end; // Generic verlet force. TGLVerletForce = class(TObject) private FOwner: TGLVerletWorld; public constructor Create(const aOwner: TGLVerletWorld); virtual; destructor Destroy; override; // Implementation should add force to force resultant for all relevant nodes procedure AddForce(const vpt: TVerletProgressTimes); virtual; abstract; // Notifies removal of a node procedure RemoveNode(const aNode: TVerletNode); virtual; abstract; property Owner: TGLVerletWorld read FOwner; end; // A verlet force that applies to two specified nodes. TGLVerletDualForce = class(TGLVerletForce) private FNodeA, FNodeB: TVerletNode; public procedure RemoveNode(const aNode: TVerletNode); override; { Reference to NodeA. } property NodeA: TVerletNode read FNodeA write FNodeA; { Reference to NodeB. } property NodeB: TVerletNode read FNodeB write FNodeB; end; // A verlet force that applies to a specified group of nodes. TVerletGroupForce = class(TGLVerletForce) private FNodes: TVerletNodeList; public constructor Create(const aOwner: TGLVerletWorld); override; destructor Destroy; override; procedure RemoveNode(const aNode: TVerletNode); override; // Nodes of the force group, referred, NOT owned. property Nodes: TVerletNodeList read FNodes; end; // A global force (applied to all verlet nodes). TGLVerletGlobalForce = class(TGLVerletForce) public procedure RemoveNode(const aNode: TVerletNode); override; procedure AddForce(const vpt: TVerletProgressTimes); override; procedure AddForceToNode(const aNode: TVerletNode); virtual; abstract; end; TGLVerletForceList = class(TList) private function GetItems(i: Integer): TGLVerletForce; procedure SetItems(i: Integer; const Value: TGLVerletForce); public property Items[i: Integer]: TGLVerletForce read GetItems write SetItems; default; end; TVCStick = class; TVFSpring = class; TVCSlider = class; TUpdateSpacePartion = (uspEveryIteration, uspEveryFrame, uspNever); TCollisionConstraintTypes = (cctEdge, cctNode); TCollisionConstraintTypesSet = set of TCollisionConstraintTypes; TGLVerletWorld = class(TObject) private FIterations: Integer; FNodes: TVerletNodeList; FConstraints: TVerletConstraintList; FForces: TGLVerletForceList; FMaxDeltaTime, FSimTime: Single; FDrag: Single; FCurrentDeltaTime: Single; FInvCurrentDeltaTime: Single; FSolidEdges: TGLVerletEdgeList; FSpacePartition: TBaseSpacePartition; FCurrentStepCount: Integer; FUpdateSpacePartion: TUpdateSpacePartion; FCollisionConstraintTypes: TCollisionConstraintTypesSet; FConstraintsWithBeforeIterations: TVerletConstraintList; FVerletNodeClass: TGLVerletNodeClass; FInertia: Boolean; FInertaPauseSteps: Integer; protected procedure AccumulateForces(const vpt: TVerletProgressTimes); virtual; procedure Verlet(const vpt: TVerletProgressTimes); virtual; procedure SatisfyConstraints(const vpt: TVerletProgressTimes); virtual; procedure DoUpdateSpacePartition; public constructor Create; virtual; destructor Destroy; override; function AddNode(const aNode: TVerletNode): Integer; procedure RemoveNode(const aNode: TVerletNode); function AddConstraint(const aConstraint: TVerletConstraint): Integer; procedure RemoveConstraint(const aConstraint: TVerletConstraint); function AddForce(const aForce: TGLVerletForce): Integer; procedure RemoveForce(const aForce: TGLVerletForce); procedure AddSolidEdge(const aNodeA, aNodeB: TVerletNode); procedure PauseInertia(const IterationSteps: Integer); function CreateOwnedNode(const Location: TAffineVector; const aRadius: Single = 0; const aWeight: Single = 1): TVerletNode; function CreateStick(const aNodeA, aNodeB: TVerletNode; const Slack: Single = 0): TVCStick; function CreateSpring(const aNodeA, aNodeB: TVerletNode; const aStrength, aDamping: Single; const aSlack: Single = 0): TVFSpring; function CreateSlider(const aNodeA, aNodeB: TVerletNode; const aSlideDirection: TAffineVector): TVCSlider; procedure Initialize; virtual; procedure CreateOctree(const OctreeMin, OctreeMax: TAffineVector; const LeafThreshold, MaxTreeDepth: Integer); function Progress(const deltaTime, newTime: Double): Integer; virtual; function FirstNode: TVerletNode; function LastNode: TVerletNode; property Drag: Single read FDrag write FDrag; property Iterations: Integer read FIterations write FIterations; property Nodes: TVerletNodeList read FNodes; property Constraints: TVerletConstraintList read FConstraints; property ConstraintsWithBeforeIterations: TVerletConstraintList read FConstraintsWithBeforeIterations; property SimTime: Single read FSimTime write FSimTime; property MaxDeltaTime: Single read FMaxDeltaTime write FMaxDeltaTime; property CurrentDeltaTime: Single read FCurrentDeltaTime; property SolidEdges: TGLVerletEdgeList read FSolidEdges write FSolidEdges; property CurrentStepCount: Integer read FCurrentStepCount; property SpacePartition: TBaseSpacePartition read FSpacePartition; property UpdateSpacePartion: TUpdateSpacePartion read FUpdateSpacePartion write FUpdateSpacePartion; property CollisionConstraintTypes: TCollisionConstraintTypesSet read FCollisionConstraintTypes write FCollisionConstraintTypes; property VerletNodeClass: TGLVerletNodeClass read FVerletNodeClass write FVerletNodeClass; property Inertia: Boolean read FInertia write FInertia; end; TVFGravity = class(TGLVerletGlobalForce) private FGravity: TAffineVector; public constructor Create(const aOwner: TGLVerletWorld); override; procedure AddForceToNode(const aNode: TVerletNode); override; property Gravity: TAffineVector read FGravity write FGravity; end; TVFAirResistance = class(TGLVerletGlobalForce) private FDragCoeff: Single; FWindDirection: TAffineVector; FWindMagnitude: Single; FWindChaos: Single; procedure SetWindDirection(const Value: TAffineVector); public constructor Create(const aOwner: TGLVerletWorld); override; procedure AddForceToNode(const aNode: TVerletNode); override; property DragCoeff: Single read FDragCoeff write FDragCoeff; property WindDirection: TAffineVector read FWindDirection write SetWindDirection; property WindMagnitude: Single read FWindMagnitude write FWindMagnitude; // Measures how chaotic the wind is, as a fraction of the wind magnitude property WindChaos: Single read FWindChaos write FWindChaos; end; TVFSpring = class(TGLVerletDualForce) private FRestLength: Single; FStrength: Single; FDamping: Single; FSlack: Single; FForceFactor: Single; protected procedure SetSlack(const Value: Single); public procedure AddForce(const vpt: TVerletProgressTimes); override; // Must be invoked after adjust node locations or strength procedure SetRestLengthToCurrent; property Strength: Single read FStrength write FStrength; property Damping: Single read FDamping write FDamping; property Slack: Single read FSlack write SetSlack; end; // Floor collision constraint TVCFloor = class(TGLVerletGlobalFrictionConstraintSP) private FBounceRatio, FFloorLevel: Single; FNormal: TAffineVector; protected procedure SetNormal(const Value: TAffineVector); public constructor Create(const aOwner: TGLVerletWorld); override; procedure PerformSpaceQuery; override; procedure SatisfyConstraintForNode(const aNode: TVerletNode; const iteration, maxIterations: Integer); override; property BounceRatio: Single read FBounceRatio write FBounceRatio; property FloorLevel: Single read FFloorLevel write FFloorLevel; property Normal: TAffineVector read FNormal write SetNormal; end; TVCHeightField = class; TVCHeightFieldOnNeedHeight = function(hfConstraint: TVCHeightField; node: TVerletNode): Single of object; // HeightField collision constraint (punctual!) TVCHeightField = class(TVCFloor) private FOnNeedHeight: TVCHeightFieldOnNeedHeight; public procedure SatisfyConstraintForNode(const aNode: TVerletNode; const iteration, maxIterations: Integer); override; property OnNeedHeight: TVCHeightFieldOnNeedHeight read FOnNeedHeight write FOnNeedHeight; end; { Stick constraint. Imposes a fixed distance between two nodes. } TVCStick = class(TGLVerletDualConstraint) private FSlack: Single; FRestLength: Single; public procedure SatisfyConstraint(const iteration, maxIterations: Integer); override; procedure SetRestLengthToCurrent; property Slack: Single read FSlack write FSlack; property RestLength: Single read FRestLength write FRestLength; end; { Rigid body constraint. Regroups several nodes in a rigid body conformation, somewhat similar to a stick but for multiple nodes. EXPERIMENTAL, DOES NOT WORK! } TVCRigidBody = class(TVerletGroupConstraint) private FNodeParams: array of TAffineVector; FNodeCoords: array of TAffineVector; FNatMatrix, FInvNatMatrix: TAffineMatrix; protected procedure ComputeBarycenter(var barycenter: TAffineVector); procedure ComputeNaturals(const barycenter: TAffineVector; var natX, natY, natZ: TAffineVector); public procedure ComputeRigidityParameters; procedure SatisfyConstraint(const iteration, maxIterations: Integer); override; end; { Slider constraint. Imposes that two nodes be aligned on a defined direction, on which they can slide freely. Note that the direction is fixed and won't rotate with the verlet assembly!. } TVCSlider = class(TGLVerletDualConstraint) private FSlideDirection: TAffineVector; FConstrained: Boolean; protected procedure SetSlideDirection(const Value: TAffineVector); public procedure SatisfyConstraint(const iteration, maxIterations: Integer); override; property SlideDirection: TAffineVector read FSlideDirection write SetSlideDirection; { Constrain NodeB to the halfplane defined by NodeA and SlideDirection. } property Constrained: Boolean read FConstrained write FConstrained; end; { Sphere collision constraint. } TVCSphere = class(TGLVerletGlobalFrictionConstraintSphere) private FRadius: Single; public function GetBSphere: TBSphere; override; procedure SatisfyConstraintForNode(const aNode: TVerletNode; const iteration, maxIterations: Integer); override; procedure SatisfyConstraintForEdge(const aEdge: TGLVerletEdge; const iteration, maxIterations: Integer); override; property Radius: Single read FRadius write FRadius; end; { Cylinder collision constraint. The cylinder is considered infinite by this constraint. } TVCCylinder = class(TGLVerletGlobalFrictionConstraint) private FAxis: TAffineVector; FRadius, FRadius2: Single; protected procedure SetRadius(const val: Single); public procedure SatisfyConstraintForNode(const aNode: TVerletNode; const iteration, maxIterations: Integer); override; { A base point on the cylinder axis. Can theoretically be anywhere, however, to reduce floating point precision issues, choose it in the area where collision detection will occur. } // property Base : TAffineVector read FBase write FBase; { Cylinder axis vector. Must be normalized. } property Axis: TAffineVector read FAxis write FAxis; { Cylinder radius. } property Radius: Single read FRadius write SetRadius; end; { Cube collision constraint. } TVCCube = class(TGLVerletGlobalFrictionConstraintBox) private FHalfSides: TAffineVector; FSides: TAffineVector; FDirection: TAffineVector; procedure SetSides(const Value: TAffineVector); public function GetAABB: TAABB; override; procedure SatisfyConstraintForNode(const aNode: TVerletNode; const iteration, maxIterations: Integer); override; // Broken and very slow! procedure SatisfyConstraintForEdge(const aEdge: TGLVerletEdge; const iteration, maxIterations: Integer); override; property Direction: TAffineVector read FDirection write FDirection; property Sides: TAffineVector read FSides write SetSides; end; // Capsule collision constraint. TVCCapsule = class(TGLVerletGlobalFrictionConstraintSphere) private FAxis: TAffineVector; FRadius, FRadius2, FLength, FLengthDiv2: Single; protected procedure SetAxis(const val: TAffineVector); procedure SetRadius(const val: Single); procedure SetLength(const val: Single); public function GetBSphere: TBSphere; override; procedure SatisfyConstraintForNode(const aNode: TVerletNode; const iteration, maxIterations: Integer); override; procedure SatisfyConstraintForEdge(const aEdge: TGLVerletEdge; const iteration, maxIterations: Integer); override; // property Base : TAffineVector read FBase write FBase; property Axis: TAffineVector read FAxis write SetAxis; property Radius: Single read FRadius write SetRadius; property Length: Single read FLength write SetLength; end; (* Specialized verlet node that can be anchored to a GLScene object. If it's anchored and has the property "NailedDown" set, it will remain in the same relative position to the GLScene object. *) TGLVerletNode = class(TVerletNode) private FRelativePosition: TAffineVector; FGLBaseSceneObject: TGLBaseSceneObject; procedure SetGLBaseSceneObject(const Value: TGLBaseSceneObject); protected procedure SetLocation(const Value: TAffineVector); override; public procedure Verlet(const vpt: TVerletProgressTimes); override; property GLBaseSceneObject: TGLBaseSceneObject read FGLBaseSceneObject write SetGLBaseSceneObject; property RelativePosition: TAffineVector read FRelativePosition write FRelativePosition; end; function CreateVCPlaneFromGLPlane(Plane: TGLPlane; VerletWorld: TGLVerletWorld; Offset: Single): TVCFloor; // ------------------------------------------------------------------ implementation // ------------------------------------------------------------------ function CreateVCPlaneFromGLPlane(Plane: TGLPlane; VerletWorld: TGLVerletWorld; Offset: Single): TVCFloor; begin result := TVCFloor.Create(VerletWorld); with result do begin Normal := VectorNormalize(Plane.Direction.AsAffineVector); Location := VectorAdd(Plane.Position.AsAffineVector, VectorScale(Normal, Offset)); end; end; // ---------------------------- { TGLVerletNode } // ---------------------------- procedure TGLVerletNode.SetGLBaseSceneObject(const Value: TGLBaseSceneObject); begin FGLBaseSceneObject := Value; if Assigned(GLBaseSceneObject) and NailedDown then FRelativePosition := AffineVectorMake(GLBaseSceneObject.AbsoluteToLocal(VectorMake(FLocation, 1))); end; procedure TGLVerletNode.SetLocation(const Value: TAffineVector); begin inherited; if Assigned(GLBaseSceneObject) and NailedDown then FRelativePosition := GLBaseSceneObject.AbsoluteToLocal(Value); end; procedure TGLVerletNode.Verlet(const vpt: TVerletProgressTimes); begin if Assigned(GLBaseSceneObject) and NailedDown then begin FLocation := GLBaseSceneObject.LocalToAbsolute(FRelativePosition); end else inherited; end; // ------------------ // ------------------ TVerletNode ------------------ // ------------------ constructor TVerletNode.CreateOwned(const aOwner: TGLVerletWorld); begin inherited CreateOwned(aOwner.SpacePartition); if Assigned(aOwner) then aOwner.AddNode(Self); FWeight := 1; FInvWeight := 1; FRadius := 0; FFriction := 1; end; destructor TVerletNode.Destroy; begin if Assigned(FOwner) then FOwner.RemoveNode(Self); inherited; end; { TODO: Improve the friction calculations Friction = - NormalForce * FrictionConstant To compute the NormalForce, which is the force acting on the normal of the collider, we can use the fact that F = m*a. m is the weight of the node, a is the acceleration (retardation) caused by the collission. Acceleration := - PenetrationDepth / Owner.FCurrentDeltaTime; The force with which the node has been "stopped" from penetration NormalForce := Weight * Acceleration; This force should be applied to stopping the movement. } procedure TVerletNode.ApplyFriction(const friction, penetrationDepth: Single; const surfaceNormal: TAffineVector); var frictionMove, move, moveNormal: TAffineVector; realFriction: Single; begin if (penetrationDepth > 0) then begin realFriction := friction * FFriction; if realFriction > 0 then begin VectorSubtract(Location, OldLocation, move); moveNormal := VectorScale(surfaceNormal, VectorDotProduct(move, surfaceNormal)); frictionMove := VectorSubtract(move, moveNormal); if penetrationDepth > Radius then ScaleVector(frictionMove, realFriction) else ScaleVector(frictionMove, realFriction * Sqrt(penetrationDepth / Radius)); VectorAdd(OldLocation, frictionMove, FOldLocation); end; end; end; procedure TVerletNode.OldApplyFriction(const friction, penetrationDepth: Single); var frictionMove, move: TAffineVector; // pd : Single; begin VectorSubtract(Location, OldLocation, move); VectorScale(move, friction * FFriction, frictionMove); // pd:=Abs(penetrationDepth); // ScaleVector(frictionMove, friction*pd); VectorAdd(OldLocation, frictionMove, FOldLocation); end; function TVerletNode.DistanceToNode(const node: TVerletNode): Single; begin Result := VectorDistance(Location, node.Location); end; function TVerletNode.GetMovement: TAffineVector; begin Result := VectorSubtract(Location, OldLocation); end; procedure TVerletNode.Initialize; begin FOldLocation := Location; end; procedure TVerletNode.SetWeight(const Value: Single); begin FWeight := Value; if Value <> 0 then FInvWeight := 1 / Value else FInvWeight := 1; end; procedure TVerletNode.Verlet(const vpt: TVerletProgressTimes); var newLocation, temp, move, accel: TAffineVector; begin if NailedDown then begin FOldLocation := Location; end else begin if Owner.Inertia then begin temp := Location; VectorSubtract(Location, OldLocation, move); ScaleVector(move, 1 - Owner.Drag); // *Sqr(deltaTime)); VectorAdd(Location, move, newLocation); VectorScale(Force, vpt.sqrDeltaTime * FInvWeight, accel); AddVector(newLocation, accel); Location := newLocation; FOldLocation := temp; end else begin newLocation := Location; VectorScale(Force, vpt.sqrDeltaTime * FInvWeight, accel); AddVector(newLocation, accel); Location := newLocation; FOldLocation := Location; end; end; end; procedure TVerletNode.AfterProgress; begin // nothing here, reserved for subclass use end; // ------------------ // ------------------ TVerletNodeList ------------------ // ------------------ function TVerletNodeList.GetItems(i: Integer): TVerletNode; begin result := Get(i); end; procedure TVerletNodeList.SetItems(i: Integer; const Value: TVerletNode); begin Put(i, Value); end; function TVerletNode.GetSpeed: TAffineVector; begin result := VectorScale(VectorSubtract(FLocation, FOldLocation), 1 / Owner.CurrentDeltaTime); end; // ------------------ // ------------------ TVerletConstraint ------------------ // ------------------ constructor TVerletConstraint.Create(const aOwner: TGLVerletWorld); begin inherited Create; if Assigned(aOwner) then aOwner.AddConstraint(Self); FEnabled := True; end; destructor TVerletConstraint.Destroy; begin if Assigned(FOwner) then FOwner.RemoveConstraint(Self); inherited; end; procedure TVerletConstraint.BeforeIterations; begin // NADA! end; // ------------------ // ------------------ TGLVerletDualConstraint ------------------ // ------------------ procedure TGLVerletDualConstraint.RemoveNode(const aNode: TVerletNode); begin if FNodeA = aNode then FNodeA := nil; if FNodeB = aNode then FNodeB := nil; if (FNodeA = nil) and (FNodeA = nil) then Free; end; // ------------------ // ------------------ TVerletGroupConstraint ------------------ // ------------------ constructor TVerletGroupConstraint.Create(const aOwner: TGLVerletWorld); begin inherited Create(aOwner); FNodes := TVerletNodeList.Create; end; destructor TVerletGroupConstraint.Destroy; begin FNodes.Free; inherited; end; procedure TVerletGroupConstraint.RemoveNode(const aNode: TVerletNode); begin FNodes.Remove(aNode); end; // ------------------ // ------------------ TGLVerletGlobalConstraint ------------------ // ------------------ procedure TGLVerletGlobalConstraint.AddKickbackForceAt(const Pos: TAffineVector; const Force: TAffineVector); var dPos: TAffineVector; begin // Sum forces AddVector(FKickbackForce, Force); // Sum torques dPos := VectorSubtract(Pos, FLocation); AddVector(FKickbackTorque, VectorCrossProduct(dPos, Force)); end; function TGLVerletGlobalConstraint.TranslateKickbackTorque(const TorqueCenter: TAffineVector): TAffineVector; begin // EM(b) = EM(a) + EF x VectorSubtract(b, a). result := VectorAdd(FKickbackTorque, VectorCrossProduct(VectorSubtract(TorqueCenter, FLocation), FKickbackForce)); end; procedure TGLVerletGlobalConstraint.BeforeIterations; begin inherited; FKickbackForce := NullVector; FKickbackTorque := NullVector; end; procedure TGLVerletGlobalConstraint.RemoveNode(const aNode: TVerletNode); begin // nothing to do here end; procedure TGLVerletGlobalConstraint.SetLocation(const Value: TAffineVector); begin FLocation := Value; end; procedure TGLVerletGlobalConstraint.SatisfyConstraint(const iteration, maxIterations: Integer); var i: Integer; node: TVerletNode; begin if cctNode in Owner.CollisionConstraintTypes then for i := 0 to Owner.Nodes.Count - 1 do begin node := TVerletNode(Owner.Nodes[i]); if not node.NailedDown then SatisfyConstraintForNode(node, iteration, maxIterations); end; // } if cctEdge in Owner.CollisionConstraintTypes then for i := 0 to Owner.SolidEdges.Count - 1 do begin SatisfyConstraintForEdge(Owner.SolidEdges[i], iteration, maxIterations); end; // } end; procedure TGLVerletGlobalConstraint.SatisfyConstraintForEdge(const aEdge: TGLVerletEdge; const iteration, maxIterations: Integer); begin // Purely virtual, but can't be abstract... end; // ------------------ // ------------------ TGLVerletGlobalFrictionConstraint ------------------ // ------------------ constructor TGLVerletGlobalFrictionConstraint.Create(const aOwner: TGLVerletWorld); begin inherited; FFrictionRatio := cDEFAULT_CONSTRAINT_FRICTION; end; // ------------------ // ------------------ TGLVerletGlobalFrictionConstraintSP ------------------ // ------------------ procedure TGLVerletGlobalFrictionConstraintSP.SatisfyConstraint(const iteration, maxIterations: Integer); var i: Integer; node: TVerletNode; edge: TGLVerletEdge; SP: TBaseSpacePartition; Leaf: TSpacePartitionLeaf; begin if Owner.SpacePartition = nil then begin inherited; Exit; end; PerformSpaceQuery; SP := Owner.SpacePartition; for i := 0 to SP.QueryResult.Count - 1 do begin Leaf := SP.QueryResult[i]; if Leaf is TVerletNode then begin if cctNode in Owner.CollisionConstraintTypes then begin node := Leaf as TVerletNode; if not node.NailedDown then SatisfyConstraintForNode(node, iteration, maxIterations); end; end else if Leaf is TGLVerletEdge then begin if cctEdge in Owner.CollisionConstraintTypes then begin edge := Leaf as TGLVerletEdge; SatisfyConstraintForEdge(edge, iteration, maxIterations); end; end else Assert(False, 'Bad objects in list!'); end; end; // ------------------ // ------------------ TVerletConstraintList ------------------ // ------------------ function TVerletConstraintList.GetItems(i: Integer): TVerletConstraint; begin result := Get(i); end; procedure TVerletConstraintList.SetItems(i: Integer; const Value: TVerletConstraint); begin Put(i, Value); end; // ------------------ // ------------------ TGLVerletForce ------------------ // ------------------ constructor TGLVerletForce.Create(const aOwner: TGLVerletWorld); begin inherited Create; if Assigned(aOwner) then aOwner.AddForce(Self); end; destructor TGLVerletForce.Destroy; begin if Assigned(FOwner) then FOwner.RemoveForce(Self); inherited; end; // ------------------ // ------------------ TVerletGroupForce ------------------ // ------------------ constructor TVerletGroupForce.Create(const aOwner: TGLVerletWorld); begin inherited Create(aOwner); FNodes := TVerletNodeList.Create; end; destructor TVerletGroupForce.Destroy; begin FNodes.Free; inherited; end; procedure TVerletGroupForce.RemoveNode(const aNode: TVerletNode); begin FNodes.Remove(aNode); end; // ------------------ // ------------------ TGLVerletGlobalForce ------------------ // ------------------ procedure TGLVerletGlobalForce.RemoveNode(const aNode: TVerletNode); begin // nothing to do here end; procedure TGLVerletGlobalForce.AddForce; var i: Integer; node: TVerletNode; begin for i := 0 to Owner.Nodes.Count - 1 do begin node := TVerletNode(Owner.Nodes.List[i]); if not node.NailedDown then AddForceToNode(node); end; end; // ------------------ // ------------------ TGLVerletDualForce ------------------ // ------------------ procedure TGLVerletDualForce.RemoveNode(const aNode: TVerletNode); begin if FNodeA = aNode then FNodeA := nil; if FNodeB = aNode then FNodeB := nil; end; // ------------------ // ------------------ TGLVerletForceList ------------------ // ------------------ function TGLVerletForceList.GetItems(i: Integer): TGLVerletForce; begin result := Get(i); end; procedure TGLVerletForceList.SetItems(i: Integer; const Value: TGLVerletForce); begin Put(i, Value); end; // ------------------ // ------------------ TGLVerletWorld ------------------ // ------------------ constructor TGLVerletWorld.Create; begin inherited; FDrag := G_DRAG; FNodes := TVerletNodeList.Create; FConstraints := TVerletConstraintList.Create; FConstraintsWithBeforeIterations := TVerletConstraintList.Create; FForces := TGLVerletForceList.Create; FMaxDeltaTime := 0.02; FIterations := 3; FSolidEdges := TGLVerletEdgeList.Create; FCurrentStepCount := 0; FUpdateSpacePartion := uspNever; FCollisionConstraintTypes := [cctNode, cctEdge]; FSpacePartition := nil; FVerletNodeClass := TVerletNode; FInertia := True; end; destructor TGLVerletWorld.Destroy; var i: Integer; begin // Delete all nodes for i := 0 to FNodes.Count - 1 do with FNodes[i] do begin FOwner := nil; Free; end; FreeAndNil(FNodes); // Delete all constraints for i := 0 to FConstraints.Count - 1 do with FConstraints[i] do begin FOwner := nil; Free; end; FreeAndNil(FConstraints); // Delete all forces for i := 0 to FForces.Count - 1 do with FForces[i] do begin FOwner := nil; Free; end; FreeAndNil(FForces); FreeAndNil(FConstraintsWithBeforeIterations); for i := 0 to FSolidEdges.Count - 1 do FSolidEdges[i].Free; FreeAndNil(FSolidEdges); FreeAndNil(FSpacePartition); inherited; end; procedure TGLVerletWorld.AccumulateForces(const vpt: TVerletProgressTimes); var i: Integer; begin // First of all, reset all forces for i := 0 to FNodes.Count - 1 do FNodes[i].FForce := NullVector; // Now, update all forces in the assembly! for i := 0 to FForces.Count - 1 do FForces[i].AddForce(vpt); end; function TGLVerletWorld.AddNode(const aNode: TVerletNode): Integer; begin if Assigned(aNode.FOwner) then aNode.Owner.FNodes.Remove(aNode); result := FNodes.Add(aNode); aNode.FOwner := Self; end; procedure TGLVerletWorld.RemoveNode(const aNode: TVerletNode); var i: Integer; begin if aNode.Owner = Self then begin FNodes.Remove(aNode); aNode.FOwner := nil; // drop refs in constraints for i := FConstraints.Count - 1 downto 0 do FConstraints[i].RemoveNode(aNode); // drop refs in forces for i := FForces.Count - 1 downto 0 do FForces[i].RemoveNode(aNode); end; end; function TGLVerletWorld.AddConstraint(const aConstraint: TVerletConstraint): Integer; begin if Assigned(aConstraint.FOwner) then aConstraint.Owner.FConstraints.Remove(aConstraint); result := FConstraints.Add(aConstraint); aConstraint.FOwner := Self; end; procedure TGLVerletWorld.RemoveConstraint(const aConstraint: TVerletConstraint); begin if aConstraint.Owner = Self then begin FConstraints.Remove(aConstraint); aConstraint.FOwner := nil; end; end; function TGLVerletWorld.AddForce(const aForce: TGLVerletForce): Integer; begin if Assigned(aForce.FOwner) then aForce.Owner.FForces.Remove(aForce); result := FForces.Add(aForce); aForce.FOwner := Self; end; procedure TGLVerletWorld.RemoveForce(const aForce: TGLVerletForce); begin if aForce.Owner = Self then begin FForces.Remove(aForce); aForce.FOwner := nil; end; end; procedure TGLVerletWorld.AddSolidEdge(const aNodeA, aNodeB: TVerletNode); var VerletEdge: TGLVerletEdge; begin VerletEdge := TGLVerletEdge.CreateEdgeOwned(aNodeA, aNodeB); SolidEdges.Add(VerletEdge); end; function TGLVerletWorld.FirstNode: TVerletNode; begin Assert(FNodes.Count > 0, 'There are no nodes in the assembly!'); result := FNodes[0]; end; function TGLVerletWorld.LastNode: TVerletNode; begin Assert(FNodes.Count > 0, 'There are no nodes in the assembly!'); result := FNodes[FNodes.Count - 1]; end; function TGLVerletWorld.CreateOwnedNode(const Location: TAffineVector; const aRadius: Single = 0; const aWeight: Single = 1): TVerletNode; begin Result:=VerletNodeClass.CreateOwned(Self); Result.Location:=Location; Result.OldLocation:=Location; Result.Weight:=aWeight; Result.Radius:=aRadius; end; function TGLVerletWorld.CreateStick(const aNodeA, aNodeB: TVerletNode; const Slack: Single = 0): TVCStick; begin Assert(aNodeA <> aNodeB, 'Can''t create stick between same node!'); Result:=TVCStick.Create(Self); Result.NodeA:=aNodeA; Result.NodeB:=aNodeB; Result.SetRestLengthToCurrent; Result.Slack := Slack; end; function TGLVerletWorld.CreateSpring(const aNodeA, aNodeB: TVerletNode; const aStrength, aDamping: Single; const aSlack: Single = 0): TVFSpring; begin Result:=TVFSpring.Create(Self); Result.NodeA:=aNodeA; Result.NodeB:=aNodeB; Result.Strength:=aStrength; Result.Damping:=aDamping; Result.Slack:=aSlack; Result.SetRestLengthToCurrent; end; function TGLVerletWorld.CreateSlider(const aNodeA, aNodeB: TVerletNode; const aSlideDirection: TAffineVector): TVCSlider; begin Result:=TVCSlider.Create(Self); Result.NodeA:=aNodeA; Result.NodeB:=aNodeB; Result.SlideDirection:=aSlideDirection; end; procedure TGLVerletWorld.Initialize; var i: Integer; begin for i := 0 to FNodes.Count - 1 do FNodes[i].Initialize; end; function TGLVerletWorld.Progress(const deltaTime, newTime: Double): Integer; var i: Integer; ticks: Integer; myDeltaTime: Single; vpt: TVerletProgressTimes; begin ticks := 0; myDeltaTime := FMaxDeltaTime; FCurrentDeltaTime := FMaxDeltaTime; FInvCurrentDeltaTime := 1 / FCurrentDeltaTime; vpt.deltaTime := myDeltaTime; vpt.sqrDeltaTime := Sqr(myDeltaTime); vpt.invSqrDeltaTime := 1 / vpt.sqrDeltaTime; while FSimTime < newTime do begin Inc(ticks); FSimTime := FSimTime + myDeltaTime; vpt.newTime := FSimTime; Verlet(vpt); AccumulateForces(vpt); SatisfyConstraints(vpt); if FInertaPauseSteps > 0 then begin dec(FInertaPauseSteps); if FInertaPauseSteps = 0 then Inertia := True; end; Break; end; result := ticks; for i := 0 to FNodes.Count - 1 do FNodes[i].AfterProgress; end; procedure TGLVerletWorld.DoUpdateSpacePartition; var i: Integer; begin if Assigned(SpacePartition) then begin for i := 0 to FSolidEdges.Count - 1 do if (FSolidEdges[i].FNodeA.FChangedOnStep = FCurrentStepCount) or (FSolidEdges[i].FNodeB.FChangedOnStep = FCurrentStepCount) then FSolidEdges[i].Changed; for i := 0 to FNodes.Count - 1 do if (FNodes[i].FChangedOnStep = FCurrentStepCount) then FNodes[i].Changed; end; end; procedure TGLVerletWorld.SatisfyConstraints(const vpt: TVerletProgressTimes); var i, j: Integer; Constraint: TVerletConstraint; begin for i := 0 to FConstraintsWithBeforeIterations.Count - 1 do begin Constraint := FConstraintsWithBeforeIterations[i]; Constraint.BeforeIterations; end; if UpdateSpacePartion = uspEveryFrame then Inc(FCurrentStepCount); for j := 0 to Iterations - 1 do begin for i := 0 to FConstraints.Count - 1 do with FConstraints[i] do if Enabled then SatisfyConstraint(j, Iterations); // } if UpdateSpacePartion = uspEveryIteration then DoUpdateSpacePartition; end; if UpdateSpacePartion = uspEveryFrame then DoUpdateSpacePartition; // } end; procedure TGLVerletWorld.Verlet(const vpt: TVerletProgressTimes); var i: Integer; begin if UpdateSpacePartion <> uspNever then Inc(FCurrentStepCount); for i := 0 to FNodes.Count - 1 do FNodes[i].Verlet(vpt); if UpdateSpacePartion <> uspNever then DoUpdateSpacePartition; end; // ------------------ // ------------------ TVFGravity ------------------ // ------------------ constructor TVFGravity.Create(const aOwner: TGLVerletWorld); begin inherited; FGravity.X := 0; FGravity.Y := -9.81; FGravity.Z := 0; end; procedure TVFGravity.AddForceToNode(const aNode: TVerletNode); begin CombineVector(aNode.FForce, Gravity, @aNode.Weight); end; // ------------------ // ------------------ TVFSpring ------------------ // ------------------ procedure TVFSpring.SetSlack(const Value: Single); begin if Value <= 0 then FSlack := 0 else FSlack := Value; end; procedure TVFSpring.AddForce; var hTerm, dTerm: Single; deltaV, Force: TAffineVector; deltaLength: Single; begin VectorSubtract(NodeA.Location, NodeB.Location, Force); deltaLength := VectorLength(Force); if deltaLength > FSlack then begin hTerm := (FRestLength - deltaLength) * FForceFactor; Force := VectorScale(Force, hTerm / deltaLength); end else Force := NullVector; if FDamping <> 0 then begin VectorSubtract(NodeA.GetMovement, NodeB.GetMovement, deltaV); dTerm := -0.25 * FDamping * vpt.invSqrDeltaTime; CombineVector(Force, deltaV, dTerm); end; AddVector(NodeA.FForce, Force); SubtractVector(NodeB.FForce, Force); end; procedure TVFSpring.SetRestLengthToCurrent; begin FRestLength := VectorDistance(NodeA.Location, NodeB.Location); FForceFactor := FStrength / FRestLength; end; // ------------------ // ------------------ TVFAirResistance ------------------ // ------------------ procedure TVFAirResistance.AddForceToNode(const aNode: TVerletNode); var s, F, FCurrentWindBurst: TAffineVector; sMag: Single; r: Single; Chaos: Single; begin s := aNode.Speed; if FWindMagnitude <> 0 then begin Chaos := FWindMagnitude * FWindChaos; FCurrentWindBurst.X := FWindDirection.X * FWindMagnitude + Chaos * (random - 0.5) * 2; FCurrentWindBurst.Y := FWindDirection.Y * FWindMagnitude + Chaos * (random - 0.5) * 2; FCurrentWindBurst.Z := FWindDirection.Z * FWindMagnitude + Chaos * (random - 0.5) * 2; s := VectorSubtract(s, FCurrentWindBurst); end; sMag := VectorLength(s); r := aNode.Radius + 1; if sMag <> 0 then begin F := VectorScale(s, -Sqr(sMag) * Sqr(r) * pi * FDragCoeff); aNode.FForce := VectorAdd(aNode.FForce, F); end; end; constructor TVFAirResistance.Create(const aOwner: TGLVerletWorld); begin inherited; FDragCoeff := 0.001; FWindDirection.X := 0; FWindDirection.Y := 0; FWindDirection.Z := 0; FWindMagnitude := 0; FWindChaos := 0; end; procedure TVFAirResistance.SetWindDirection(const Value: TAffineVector); begin FWindDirection := VectorNormalize(Value); end; // ------------------ // ------------------ TVCFloor ------------------ // ------------------ constructor TVCFloor.Create(const aOwner: TGLVerletWorld); begin inherited; MakeVector(FNormal, 0, 1, 0); MakeVector(FLocation, 0, 0, 0); end; procedure TVCFloor.PerformSpaceQuery; begin Owner.SpacePartition.QueryPlane(FLocation, FNormal); end; procedure TVCFloor.SatisfyConstraintForNode(const aNode: TVerletNode; const iteration, maxIterations: Integer); var penetrationDepth: Single; currentPenetrationDepth: Single; d: TAffineVector; correction: TAffineVector; begin currentPenetrationDepth := -PointPlaneDistance(aNode.Location, FLocation, FNormal) + aNode.Radius + FFloorLevel; // Record how far down the node goes penetrationDepth := currentPenetrationDepth; // Correct the node location if currentPenetrationDepth > 0 then begin correction := VectorScale(FNormal, currentPenetrationDepth); if BounceRatio > 0 then begin d := VectorSubtract(aNode.FLocation, aNode.FOldLocation); if FrictionRatio > 0 then aNode.ApplyFriction(FrictionRatio, penetrationDepth, FNormal); AddVector(aNode.FLocation, correction); aNode.FOldLocation := VectorAdd(aNode.FLocation, VectorScale(d, BounceRatio)); end else begin AddVector(aNode.FLocation, correction); if FrictionRatio > 0 then aNode.ApplyFriction(FrictionRatio, penetrationDepth, FNormal); aNode.FChangedOnStep := Owner.CurrentStepCount; end; end; end; procedure TVCFloor.SetNormal(const Value: TAffineVector); begin FNormal := Value; NormalizeVector(FNormal); end; // ------------------ // ------------------ TVCHeightField ------------------ // ------------------ procedure TVCHeightField.SatisfyConstraintForNode(const aNode: TVerletNode; const iteration, maxIterations: Integer); var penetrationDepth: Single; currentPenetrationDepth: Single; d: TAffineVector; correction: TAffineVector; begin currentPenetrationDepth := -PointPlaneDistance(aNode.Location, FLocation, FNormal) + aNode.Radius; if Assigned(FOnNeedHeight) then currentPenetrationDepth := currentPenetrationDepth + FOnNeedHeight(Self, aNode); // Record how far down the node goes penetrationDepth := currentPenetrationDepth; // Correct the node location if currentPenetrationDepth > 0 then begin correction := VectorScale(FNormal, currentPenetrationDepth); if BounceRatio > 0 then begin d := VectorSubtract(aNode.FLocation, aNode.FOldLocation); if FrictionRatio > 0 then aNode.ApplyFriction(FrictionRatio, penetrationDepth, FNormal); AddVector(aNode.FLocation, correction); aNode.FOldLocation := VectorAdd(aNode.FLocation, VectorScale(d, BounceRatio)); end else begin AddVector(aNode.FLocation, correction); if FrictionRatio > 0 then aNode.ApplyFriction(FrictionRatio, penetrationDepth, FNormal); aNode.FChangedOnStep := Owner.CurrentStepCount; end; end; end; // ------------------ // ------------------ TVCStick ------------------ // ------------------ procedure TVCStick.SatisfyConstraint(const iteration, maxIterations: Integer); var delta: TAffineVector; F, r: Single; deltaLength, diff: Single; const cDefaultDelta: TAffineVector = (X: 0.01; Y: 0; Z: 0); begin Assert((NodeA <> NodeB), 'The nodes are identical - that causes division by zero!'); VectorSubtract(NodeB.Location, NodeA.Location, delta); deltaLength := VectorLength(delta); // Avoid div by zero! if deltaLength < 1E-3 then begin delta := cDefaultDelta; deltaLength := 0.01; end; diff := (deltaLength - RestLength) / deltaLength; if Abs(diff) > Slack then begin r := 1 / (NodeA.InvWeight + NodeB.InvWeight); if diff < 0 then diff := (diff + Slack) * r else diff := (diff - Slack) * r; // Take into acount the different weights of the nodes! if not NodeA.NailedDown then begin F := diff * NodeA.InvWeight; CombineVector(NodeA.FLocation, delta, F); NodeA.FChangedOnStep := Owner.CurrentStepCount; end; if not NodeB.NailedDown then begin F := -diff * NodeB.InvWeight; CombineVector(NodeB.FLocation, delta, F); NodeB.FChangedOnStep := Owner.CurrentStepCount; end; end; end; procedure TVCStick.SetRestLengthToCurrent; begin FRestLength := VectorDistance(NodeA.Location, NodeB.Location); end; // ------------------ // ------------------ TVCRigidBody ------------------ // ------------------ procedure TVCRigidBody.ComputeBarycenter(var barycenter: TAffineVector); var i: Integer; totWeight: Single; begin // first we compute the barycenter totWeight := 0; barycenter := NullVector; for i := 0 to Nodes.Count - 1 do with Nodes[i] do begin CombineVector(barycenter, Location, @Weight); totWeight := totWeight + Weight; end; if totWeight > 0 then ScaleVector(barycenter, 1 / totWeight); end; procedure TVCRigidBody.ComputeNaturals(const barycenter: TAffineVector; var natX, natY, natZ: TAffineVector); var i: Integer; delta: TAffineVector; begin natX := NullVector; natY := NullVector; natZ := NullVector; for i := 0 to Nodes.Count - 1 do begin delta := VectorSubtract(Nodes[i].Location, barycenter); CombineVector(natX, delta, FNodeParams[i].X); CombineVector(natY, delta, FNodeParams[i].Y); CombineVector(natZ, delta, FNodeParams[i].Z); end; end; procedure TVCRigidBody.ComputeRigidityParameters; var i: Integer; barycenter: TAffineVector; d: Single; begin // first we compute the barycenter ComputeBarycenter(barycenter); // next the parameters SetLength(FNodeParams, Nodes.Count); SetLength(FNodeCoords, Nodes.Count); for i := 0 to Nodes.Count - 1 do begin FNodeCoords[i] := VectorSubtract(Nodes[i].Location, barycenter); d := Nodes[i].Weight / VectorLength(FNodeCoords[i]); FNodeParams[i].X := FNodeCoords[i].X * d; FNodeParams[i].Y := FNodeCoords[i].Y * d; FNodeParams[i].Z := FNodeCoords[i].Z * d; end; ComputeNaturals(barycenter, FNatMatrix.X, FNatMatrix.Y, FNatMatrix.Z); FNatMatrix.Z:=VectorCrossProduct(FNatMatrix.X, FNatMatrix.Y); FNatMatrix.Y:=VectorCrossProduct(FNatMatrix.Z, FNatMatrix.X); NormalizeVector(FNatMatrix.X); NormalizeVector(FNatMatrix.Y); NormalizeVector(FNatMatrix.Z); FInvNatMatrix := FNatMatrix; // TransposeMatrix(FInvNatMatrix); InvertMatrix(FInvNatMatrix); end; procedure TVCRigidBody.SatisfyConstraint(const iteration, maxIterations: Integer); var i: Integer; barycenter, delta: TAffineVector; nrjBase, nrjAdjust: TAffineVector; natural: array [0 .. 2] of TAffineVector; deltas: array of TAffineVector; begin Assert(Nodes.Count = Length(FNodeParams), 'You forgot to call ComputeRigidityParameters!'); // compute the barycenter ComputeBarycenter(barycenter); // compute the natural axises ComputeNaturals(barycenter, natural[0], natural[1], natural[2]); natural[2] := VectorCrossProduct(natural[0], natural[1]); natural[1] := VectorCrossProduct(natural[2], natural[0]); for i := 0 to 2 do NormalizeVector(natural[i]); natural[0] := VectorTransform(natural[0], FInvNatMatrix); natural[1] := VectorTransform(natural[1], FInvNatMatrix); natural[2] := VectorTransform(natural[2], FInvNatMatrix); // make the natural axises orthonormal, by picking the longest two { for i:=0 to 2 do vectNorm[i]:=VectorNorm(natural[i]); if (vectNorm[0] NodeB), 'The nodes are identical - that causes division by zero!'); // project B in the plane defined by A and SlideDirection projB := VectorSubtract(NodeB.Location, NodeA.Location); F := VectorDotProduct(projB, SlideDirection); projB := VectorCombine(NodeB.Location, SlideDirection, 1, -F); if Constrained and (F < 0) then NodeB.Location := projB; VectorSubtract(projB, NodeA.Location, delta); // Take into acount the different weights of the nodes! r := 1 / (NodeA.InvWeight + NodeB.InvWeight); if not NodeA.NailedDown then begin F := r * NodeA.InvWeight; CombineVector(NodeA.FLocation, delta, F); NodeA.FChangedOnStep := Owner.CurrentStepCount; end; if not NodeB.NailedDown then begin F := -r * NodeB.InvWeight; CombineVector(NodeB.FLocation, delta, F); NodeB.FChangedOnStep := Owner.CurrentStepCount; end; end; // ------------------ // ------------------ TVCSphere ------------------ // ------------------ function TVCSphere.GetBSphere: TBSphere; begin result.Center := FLocation; result.Radius := FRadius; end; procedure TVCSphere.SatisfyConstraintForEdge(const aEdge: TGLVerletEdge; const iteration, maxIterations: Integer); var closestPoint, move, delta, contactNormal: TAffineVector; deltaLength, diff: Single; begin // If the edge penetrates the sphere, try pushing the nodes until it no // longer does closestPoint := PointSegmentClosestPoint(FLocation, aEdge.NodeA.FLocation, aEdge.NodeB.FLocation); // Find the distance between the two VectorSubtract(closestPoint, Location, delta); deltaLength := VectorLength(delta); if deltaLength < Radius then begin if deltaLength > 0 then begin contactNormal := VectorScale(delta, 1 / deltaLength); aEdge.NodeA.ApplyFriction(FFrictionRatio, Radius - Abs(deltaLength), contactNormal); aEdge.NodeB.ApplyFriction(FFrictionRatio, Radius - Abs(deltaLength), contactNormal); end; // Move it outside the sphere! diff := (Radius - deltaLength) / deltaLength; VectorScale(delta, diff, move); AddVector(aEdge.NodeA.FLocation, move); AddVector(aEdge.NodeB.FLocation, move); // Add the force to the kickback // F = a * m // a = move / deltatime AddKickbackForceAt(FLocation, VectorScale(move, -(aEdge.NodeA.FWeight + aEdge.NodeB.FWeight) * Owner.FInvCurrentDeltaTime)); aEdge.NodeA.FChangedOnStep := Owner.CurrentStepCount; aEdge.NodeB.FChangedOnStep := Owner.CurrentStepCount; end; end; procedure TVCSphere.SatisfyConstraintForNode(const aNode: TVerletNode; const iteration, maxIterations: Integer); var delta, move, contactNormal: TAffineVector; deltaLength, diff: Single; begin // Find the distance between the two VectorSubtract(aNode.Location, Location, delta); // Is it inside the sphere? deltaLength := VectorLength(delta) - aNode.Radius; if Abs(deltaLength) < Radius then begin if deltaLength > 0 then begin contactNormal := VectorScale(delta, 1 / deltaLength); aNode.ApplyFriction(FFrictionRatio, Radius - Abs(deltaLength), contactNormal); end else // Slow it down - this part should not be fired aNode.OldApplyFriction(FFrictionRatio, Radius - Abs(deltaLength)); // Move it outside the sphere! diff := (Radius - deltaLength) / deltaLength; VectorScale(delta, diff, move); AddVector(aNode.FLocation, move); aNode.FChangedOnStep := Owner.CurrentStepCount; // Add the force to the kickback // F = a * m // a = move / deltatime AddKickbackForceAt(FLocation, VectorScale(move, -aNode.FWeight * Owner.FInvCurrentDeltaTime)); end; end; // ------------------ // ------------------ TVCCylinder ------------------ // ------------------ procedure TVCCylinder.SetRadius(const val: Single); begin FRadius := val; FRadius2 := Sqr(val); end; procedure TVCCylinder.SatisfyConstraintForNode(const aNode: TVerletNode; const iteration, maxIterations: Integer); var proj, newLocation, move: TAffineVector; F, dist2, penetrationDepth: Single; begin // Compute projection of node position on the axis F := PointProject(aNode.Location, FLocation, FAxis); proj := VectorCombine(FLocation, FAxis, 1, F); // Sqr distance dist2 := VectorDistance2(proj, aNode.Location); if dist2 < FRadius2 then begin // move out of the cylinder VectorLerp(proj, aNode.Location, FRadius * RSqrt(dist2), newLocation); move := VectorSubtract(aNode.FLocation, newLocation); penetrationDepth := VectorLength(move); aNode.ApplyFriction(FFrictionRatio, penetrationDepth, VectorScale(move, 1 / penetrationDepth)); aNode.FLocation := newLocation; aNode.FChangedOnStep := Owner.CurrentStepCount; end; end; // ------------------ // ------------------ TVCCube ------------------ // ------------------ function TVCCube.GetAABB: TAABB; begin VectorAdd(FLocation, FHalfSides, result.max); VectorSubtract(FLocation, FHalfSides, result.min); end; // BROKEN AND VERY SLOW! procedure TVCCube.SatisfyConstraintForEdge(const aEdge: TGLVerletEdge; const iteration, maxIterations: Integer); var Corners: array [0 .. 7] of TAffineVector; EdgeRelative: array [0 .. 1] of TAffineVector; shortestMove { , contactNormal } : TAffineVector; shortestDeltaLength: Single; procedure AddCorner(CornerID: Integer; X, Y, Z: Single); begin X := (X - 0.5) * 2; Y := (Y - 0.5) * 2; Z := (Z - 0.5) * 2; MakeVector(Corners[CornerID], FHalfSides.X * X, FHalfSides.Y * Y, FHalfSides.Z * Z); AddVector(Corners[CornerID], FLocation); end; procedure TryEdge(Corner0, Corner1: Integer); var CubeEdgeClosest, aEdgeClosest: TAffineVector; CenteraEdge, move: TAffineVector; deltaLength: Single; begin SegmentSegmentClosestPoint( Corners[Corner0], Corners[Corner1], aEdge.NodeA.FLocation, aEdge.NodeB.FLocation, CubeEdgeClosest, aEdgeClosest); CenteraEdge := VectorSubtract(aEdgeClosest, FLocation); if (Abs(CenteraEdge.X) < FHalfSides.X) and (Abs(CenteraEdge.Y) < FHalfSides.Y) and (Abs(CenteraEdge.Z) < FHalfSides.Z) then begin // The distance to move the edge is the difference between CenterCubeEdge and // CenteraEdge move := VectorSubtract(CubeEdgeClosest, aEdgeClosest); deltaLength := VectorLength(move); if (deltaLength > 0) and (deltaLength < shortestDeltaLength) then begin shortestDeltaLength := deltaLength; shortestMove := move; end; end; end; begin // DISABLED! Exit; // Early out test EdgeRelative[0] := VectorSubtract(aEdge.FNodeA.FLocation, FLocation); EdgeRelative[1] := VectorSubtract(aEdge.FNodeB.FLocation, FLocation); // If both edges are on the same side of _any_ box side, the edge can't // cut the box if ((EdgeRelative[0].X > FHalfSides.X) and (EdgeRelative[1].X > FHalfSides.X)) or ((EdgeRelative[0].X < -FHalfSides.X) and (EdgeRelative[1].X < -FHalfSides.X)) or ((EdgeRelative[0].Y > FHalfSides.Y) and (EdgeRelative[1].Y > FHalfSides.Y)) or ((EdgeRelative[0].Y < -FHalfSides.Y) and (EdgeRelative[1].Y < -FHalfSides.Y)) or ((EdgeRelative[0].Z > FHalfSides.Z) and (EdgeRelative[1].Z > FHalfSides.Z)) or ((EdgeRelative[0].Z < -FHalfSides.Z) and (EdgeRelative[1].Z < -FHalfSides.Z)) then begin Exit; end; // For each cube edge: // find closest positions between CubeEdge and aEdge // if aEdgeClosestPosition within cube then // move nodes until closest position is outside cube // exit AddCorner(0, 0, 0, 0); AddCorner(1, 1, 0, 0); AddCorner(2, 1, 1, 0); AddCorner(3, 0, 1, 0); AddCorner(4, 0, 0, 1); AddCorner(5, 1, 0, 1); AddCorner(6, 1, 1, 1); AddCorner(7, 0, 1, 1); shortestDeltaLength := 10E30; TryEdge(0, 1); TryEdge(1, 2); TryEdge(2, 3); TryEdge(3, 0); TryEdge(4, 5); TryEdge(5, 6); TryEdge(6, 7); TryEdge(7, 4); TryEdge(0, 3); TryEdge(1, 5); TryEdge(2, 6); TryEdge(3, 7); if shortestDeltaLength < 10E8 then begin // contactNormal := VectorScale(shortestMove, 1/shortestDeltaLength); { aEdge.NodeA.ApplyFriction(FFrictionRatio, shortestDeltaLength, contactNormal); aEdge.NodeB.ApplyFriction(FFrictionRatio, shortestDeltaLength, contactNormal);// } AddVector(aEdge.NodeA.FLocation, shortestMove); AddVector(aEdge.NodeB.FLocation, shortestMove); aEdge.NodeA.Changed; aEdge.NodeB.Changed; aEdge.NodeA.FChangedOnStep := Owner.CurrentStepCount; aEdge.NodeB.FChangedOnStep := Owner.CurrentStepCount; end; end; procedure TVCCube.SatisfyConstraintForNode(const aNode: TVerletNode; const iteration, maxIterations: Integer); var p, absP, contactNormal: TAffineVector; dp: Single; smallestSide: Integer; begin // TODO: Direction of Cube should be used to rotate the nodes location, as it // stands, the cube can only face in one direction. p := VectorSubtract(aNode.FLocation, FLocation); absP.X := FHalfSides.X - Abs(p.X); absP.Y := FHalfSides.Y - Abs(p.Y); absP.Z := FHalfSides.Z - Abs(p.Z); if (PInteger(@absP.X)^ <= 0) or (PInteger(@absP.Y)^ <= 0) or (PInteger(@absP.Z)^ <= 0) then Exit; if absP.X < absP.Y then if absP.X < absP.Z then smallestSide := 0 else smallestSide := 2 else if absP.Y < absP.Z then smallestSide := 1 else smallestSide := 2; contactNormal := NullVector; // Only move along the "shortest" axis if PInteger(@p.V[smallestSide])^ >= 0 then begin dp := absP.V[smallestSide]; contactNormal.V[smallestSide] := 1; aNode.ApplyFriction(FFrictionRatio, dp, contactNormal); aNode.FLocation.V[smallestSide] := aNode.FLocation.V[smallestSide] + dp; end else begin dp := absP.V[smallestSide]; contactNormal.V[smallestSide] := -1; aNode.ApplyFriction(FFrictionRatio, dp, contactNormal); aNode.FLocation.V[smallestSide] := aNode.FLocation.V[smallestSide] - dp; end; aNode.FChangedOnStep := Owner.CurrentStepCount; end; procedure TVCCube.SetSides(const Value: TAffineVector); begin FSides := Value; FHalfSides := VectorScale(Sides, 0.5); UpdateCachedAABB; end; // ------------------ // ------------------ TVCCapsule ------------------ // ------------------ procedure TVCCapsule.SetAxis(const val: TAffineVector); begin FAxis := VectorNormalize(val); UpdateCachedBSphere; end; procedure TVCCapsule.SetLength(const val: Single); begin FLength := val; FLengthDiv2 := val * 0.5; UpdateCachedBSphere; end; procedure TVCCapsule.SetRadius(const val: Single); begin FRadius := val; FRadius2 := Sqr(val); UpdateCachedBSphere; end; function TVCCapsule.GetBSphere: TBSphere; begin result.Center := FLocation; result.Radius := Length + Radius; end; procedure TVCCapsule.SatisfyConstraintForNode(const aNode: TVerletNode; const iteration, maxIterations: Integer); var p, n2, penetrationDepth: Single; closest, v: TAffineVector; newLocation, move: TAffineVector; begin // Find the closest point to location on the capsule axis p := ClampValue(PointProject(aNode.Location, FLocation, FAxis), -FLengthDiv2, FLengthDiv2); closest := VectorCombine(FLocation, FAxis, 1, p); // vector from closest to location VectorSubtract(aNode.Location, closest, v); // should it be altered? n2 := VectorNorm(v); if n2 < FRadius2 then begin newLocation := VectorCombine(closest, v, 1, Sqrt(FRadius2 / n2)); // Do friction calculations move := VectorSubtract(newLocation, aNode.FLocation); penetrationDepth := VectorLength(move); // aNode.OldApplyFriction(FFrictionRatio, penetrationDepth); aNode.ApplyFriction(FFrictionRatio, penetrationDepth, VectorScale(move, 1 / penetrationDepth)); aNode.FLocation := newLocation; aNode.FChangedOnStep := Owner.CurrentStepCount; AddKickbackForceAt(FLocation, VectorScale(move, -aNode.FWeight * Owner.FInvCurrentDeltaTime)); end; end; procedure TVCCapsule.SatisfyConstraintForEdge(const aEdge: TGLVerletEdge; const iteration, maxIterations: Integer); var sphereLocation, closestPoint, dummy, delta, move, contactNormal: TAffineVector; Ax0, Ax1: TAffineVector; deltaLength, diff, penetrationDepth: Single; begin VectorScale(FAxis, FLengthDiv2, Ax0); AddVector(Ax0, FLocation); VectorScale(FAxis, -FLengthDiv2, Ax1); AddVector(Ax1, FLocation); SegmentSegmentClosestPoint(aEdge.NodeA.FLocation, aEdge.NodeB.FLocation, Ax0, Ax1, dummy, sphereLocation); // If the edge penetrates the sphere, try pushing the nodes until it no // longer does closestPoint := PointSegmentClosestPoint(sphereLocation, aEdge.NodeA.FLocation, aEdge.NodeB.FLocation); // Find the distance between the two VectorSubtract(closestPoint, sphereLocation, delta); deltaLength := VectorLength(delta); if deltaLength < Radius then begin // Move it outside the sphere! diff := (Radius - deltaLength) / deltaLength; VectorScale(delta, diff, move); penetrationDepth := VectorLength(move); contactNormal := VectorScale(move, 1 / penetrationDepth); aEdge.NodeA.ApplyFriction(FFrictionRatio, penetrationDepth, contactNormal); aEdge.NodeB.ApplyFriction(FFrictionRatio, penetrationDepth, contactNormal); AddVector(aEdge.NodeA.FLocation, move); AddVector(aEdge.NodeB.FLocation, move); aEdge.NodeA.FChangedOnStep := Owner.CurrentStepCount; aEdge.NodeB.FChangedOnStep := Owner.CurrentStepCount; AddKickbackForceAt(FLocation, VectorScale(move, -(aEdge.NodeA.FWeight + aEdge.NodeB.FWeight) * Owner.FInvCurrentDeltaTime)); end; end; // ------------------ // ------------------ TGLVerletEdge ------------------ // ------------------ constructor TGLVerletEdge.CreateEdgeOwned(const aNodeA, aNodeB: TVerletNode); begin FNodeA := aNodeA; FNodeB := aNodeB; inherited CreateOwned(aNodeA.Owner.SpacePartition); end; procedure TGLVerletEdge.UpdateCachedAABBAndBSphere; begin FCachedAABB.min := FNodeA.FLocation; FCachedAABB.max := FNodeA.FLocation; AABBInclude(FCachedAABB, FNodeB.FLocation); AABBToBSphere(FCachedAABB, FCachedBSphere); end; // ------------------ // ------------------ TGLVerletEdgeList ------------------ // ------------------ function TGLVerletEdgeList.GetItems(i: Integer): TGLVerletEdge; begin result := Get(i); end; procedure TGLVerletEdgeList.SetItems(i: Integer; const Value: TGLVerletEdge); begin Put(i, Value); end; procedure TVerletNode.UpdateCachedAABBAndBSphere; begin FCachedAABB.min := FLocation; FCachedAABB.max := FLocation; FCachedBSphere.Center := FLocation; FCachedBSphere.Radius := 0; end; procedure TVerletNode.SetLocation(const Value: TAffineVector); begin FLocation := Value; FChangedOnStep := Owner.CurrentStepCount; end; procedure TGLVerletWorld.CreateOctree(const OctreeMin, OctreeMax: TAffineVector; const LeafThreshold, MaxTreeDepth: Integer); var Octree: TOctreeSpacePartition; begin Assert(FNodes.Count = 0, 'You can only create an octree while the world is empty!'); FreeAndNil(FSpacePartition); Octree := TOctreeSpacePartition.Create; Octree.SetSize(OctreeMin, OctreeMax); Octree.MaxTreeDepth := MaxTreeDepth; Octree.LeafThreshold := LeafThreshold; Octree.CullingMode := cmGrossCulling; FSpacePartition := Octree; if FUpdateSpacePartion = uspNever then FUpdateSpacePartion := uspEveryFrame; end; procedure TGLVerletWorld.PauseInertia(const IterationSteps: Integer); begin FInertaPauseSteps := IterationSteps + 1; Inertia := False; end; { TGLVerletGlobalFrictionConstraintBox } procedure TGLVerletGlobalFrictionConstraintBox.PerformSpaceQuery; begin Owner.SpacePartition.QueryAABB(FCachedAABB); end; procedure TGLVerletGlobalFrictionConstraintBox.SetLocation(const Value: TAffineVector); begin inherited; UpdateCachedAABB; end; procedure TGLVerletGlobalFrictionConstraintBox.UpdateCachedAABB; begin FCachedAABB := GetAABB; end; { TGLVerletGlobalFrictionConstraintSphere } procedure TGLVerletGlobalFrictionConstraintSphere.PerformSpaceQuery; begin Owner.SpacePartition.QueryBSphere(FCachedBSphere); end; procedure TGLVerletGlobalFrictionConstraintSphere.SetLocation(const Value: TAffineVector); begin inherited; UpdateCachedBSphere; end; procedure TGLVerletGlobalFrictionConstraintSphere.UpdateCachedBSphere; begin FCachedBSphere := GetBSphere; end; constructor TGLVerletGlobalConstraint.Create(const aOwner: TGLVerletWorld); begin inherited; if Assigned(aOwner) then aOwner.ConstraintsWithBeforeIterations.Add(Self); end; destructor TGLVerletGlobalConstraint.Destroy; begin if Assigned(Owner) then Owner.ConstraintsWithBeforeIterations.Remove(Self); inherited; end; end.