From 719b38682262ccd0df5d8183b72eec5c3455fd85 Mon Sep 17 00:00:00 2001 From: iegod Date: Wed, 3 Dec 2025 10:21:36 -0500 Subject: [PATCH] Big push to implement 10mp fluid simulations. --- elements/fluidsim10.go | 77 +++++ elements/fluidsimd.go | 443 +++++++++++++++++++++++++ elements/mappedentity.go | 18 + elements/mappedentitybase.go | 51 +++ fluid/field.go | 74 +++++ fluid/fieldvector.go | 36 ++ fluid/fluid.go | 626 +++++++++++++++++++++++++++++++++++ fluid/particle.go | 6 + game/game.go | 509 ++++++++++------------------ gamedata/vector.go | 22 ++ licenses/TenMinutePhysics | 11 + 11 files changed, 1536 insertions(+), 337 deletions(-) create mode 100644 elements/fluidsim10.go create mode 100644 elements/fluidsimd.go create mode 100644 elements/mappedentity.go create mode 100644 elements/mappedentitybase.go create mode 100644 fluid/field.go create mode 100644 fluid/fieldvector.go create mode 100644 fluid/fluid.go create mode 100644 fluid/particle.go create mode 100644 licenses/TenMinutePhysics diff --git a/elements/fluidsim10.go b/elements/fluidsim10.go new file mode 100644 index 0000000..4d62732 --- /dev/null +++ b/elements/fluidsim10.go @@ -0,0 +1,77 @@ +package elements + +import ( + "fluids/fluid" + "image/color" + + "github.com/hajimehoshi/ebiten/v2" + "github.com/hajimehoshi/ebiten/v2/vector" +) + +const ( + FS10PixelWidth = 200 + FS10PixelHeight = 100 + FS10FluidWidth = 3.0 //meters + FS10FluidHeight = 1.0 //meters + FS10Resolution = 20 //need to workshop this +) + +type FluidSim10 struct { + MappedEntityBase + fluid *fluid.Fluid + angle float64 + particlebuff *ebiten.Image +} + +func NewFluidSim10() *FluidSim10 { + fsim := &FluidSim10{ + fluid: fluid.NewFluid(fluid.FieldVector{X: FS10FluidWidth, Y: FS10FluidHeight}, FS10FluidHeight/FS10Resolution), + particlebuff: ebiten.NewImage(1, 1), + } + fsim.Initialize() + return fsim +} + +func (f *FluidSim10) Initialize() { + f.Sprite = ebiten.NewImage(FS10PixelWidth, FS10PixelHeight) + f.particlebuff.Fill(color.White) + f.fluid.Initialize() +} + +func (f *FluidSim10) SetAngle(angle float64) { + f.angle = angle + f.fluid.SetAngle(float32(angle)) +} + +func (f *FluidSim10) GetAngle() float64 { + return f.angle +} + +func (f *FluidSim10) Draw() { + f.Sprite.Clear() + + vector.StrokeRect(f.Sprite, 0, 0, FS10PixelWidth, FS10PixelHeight, 2, color.White, true) + + for i := range f.fluid.Particles { + //for each particle, compute its relative position based on its + //position within the fluid field + p := &f.fluid.Particles[i] + percentx := p.Position.X / f.fluid.Field.Dimensions.X + percenty := p.Position.Y / f.fluid.Field.Dimensions.Y + ox := float64(percentx * FS10PixelWidth) + oy := float64(percenty * FS10PixelHeight) + + op := &ebiten.DrawImageOptions{} + op.GeoM.Translate(ox, oy) + f.Sprite.DrawImage(f.particlebuff, op) + } +} + +func (f *FluidSim10) Update() { + + if f.paused { + return + } + + f.fluid.Step() +} diff --git a/elements/fluidsimd.go b/elements/fluidsimd.go new file mode 100644 index 0000000..1252117 --- /dev/null +++ b/elements/fluidsimd.go @@ -0,0 +1,443 @@ +package elements + +import ( + "fluids/gamedata" + "fluids/quadtree" + "fmt" + "image/color" + "math" + "math/rand/v2" + + "github.com/hajimehoshi/ebiten/v2" + "github.com/hajimehoshi/ebiten/v2/vector" +) + +const ( + FSDWidth = 640 + FSDHeight = 360 + FSDParticleCount = 2000 + FSDGravity = 2 + FSDParticleRadius = 2.5 + FSDDamping = .7 + FSDDeltaTimeStep = 0.5 + FSDInfluenceRadius = 30 +) + +type FluidSimD struct { + MappedEntityBase + particles []*Particle + cycle int + particlebox *gamedata.Vector + particlebuff *ebiten.Image + quadtree *quadtree.Quadtree + collisionquad quadtree.Quadrant + paused bool + renderquads bool + resolvecollisions bool + resolvers []func(particle *Particle) + resolveridx int + angle float64 +} + +func NewFluidSimD() *FluidSimD { + fsd := &FluidSimD{ + particlebuff: ebiten.NewImage(FSDParticleRadius*2, FSDParticleRadius*2), + paused: false, + renderquads: false, + resolvecollisions: false, + resolveridx: 0, + angle: 0, + } + fsd.dimensions = gamedata.Vector{X: FSDWidth, Y: FSDHeight} + + //prepare root quadtree and subsequent collision search quadrant + quad := quadtree.Quadrant{ + Position: gamedata.Vector{X: FSDWidth / 2, Y: FSDHeight / 2}, + Dimensions: gamedata.Vector{X: FSDWidth, Y: FSDHeight}, + } + fsd.quadtree = quadtree.New(quad, 0) + + fsd.collisionquad = quadtree.Quadrant{ + Position: gamedata.Vector{ + X: 0, + Y: 0, + }, + Dimensions: gamedata.Vector{ + X: FSDParticleRadius * 2, + Y: FSDParticleRadius * 2, + }, + } + + //add all resolvers to strategy list + fsd.resolvers = append(fsd.resolvers, fsd.ResolveCollisionsA) + fsd.resolvers = append(fsd.resolvers, fsd.ResolveCollisionsB) + fsd.resolvers = append(fsd.resolvers, fsd.ResolveCollisionsC) + + //initialize particles, set bounding box, prepare image buffer + fsd.Sprite = ebiten.NewImage(FSDWidth, FSDHeight) + fsd.particlebox = &gamedata.Vector{ + X: FSDWidth - 50, + Y: FSDHeight - 50, + } + fsd.InitializeParticles() + + vector.FillCircle(fsd.particlebuff, FSDParticleRadius, FSDParticleRadius, FSDParticleRadius, color.White, true) + + return fsd + +} + +func (f *FluidSimD) Draw() { + f.Sprite.Clear() + + f.RenderParticles() + f.RenderBox() + + if f.renderquads { + f.RenderQuadrants() + } +} + +func (f *FluidSimD) Update() { + + if f.paused { + return + } + + f.RebuildQuadtree() + f.UpdateParticles() + + f.cycle++ +} + +func (f *FluidSimD) UpdateParticles() { + + dt := FSDDeltaTimeStep + + mx, my := ebiten.CursorPosition() + mpos := gamedata.Vector{X: float64(mx), Y: float64(my)} + maxdeflect := 40 * dt * FSDGravity + + gravity := gamedata.Vector{ + X: FSDGravity * dt * math.Sin(f.angle), + Y: FSDGravity * dt * math.Cos(f.angle), + } + + for _, particle := range f.particles { + //particle.Velocity.Y += FSDGravity * dt + + particle.Velocity = particle.Velocity.Add(gravity) + + //if inpututil.IsMouseButtonJustPressed(ebiten.MouseButtonLeft) { + + if ebiten.IsMouseButtonPressed(ebiten.MouseButtonLeft) { + delta := gamedata.Vector{ + X: mpos.X - particle.Position.X, + Y: mpos.Y - particle.Position.Y, + } + + dist := math.Sqrt(delta.X*delta.X + delta.Y*delta.Y) + theta := math.Atan2(delta.Y, delta.X) + + if dist < FSDInfluenceRadius { + + dx := dist * math.Cos(theta) + dy := dist * math.Sin(theta) + + if dx != 0 { + gainx := (-1./FSDInfluenceRadius)*math.Abs(dx) + 1. + particle.Velocity.X += 20 * gainx * -1 * math.Copysign(1, dx) + } + + if dy != 0 { + gainy := (-1./FSDInfluenceRadius)*math.Abs(dy) + 1. + particle.Velocity.Y += maxdeflect * gainy * -1 * math.Copysign(1, dy) + } + } + + } + + particle.Position.X += particle.Velocity.X * dt + particle.Position.Y += particle.Velocity.Y * dt + + if f.resolvecollisions { + f.resolvers[f.resolveridx](particle) + } + } + + for _, p := range f.particles { + f.BoundParticle(p) + } + +} + +func (f *FluidSimD) InitializeParticles() { + + f.particles = f.particles[:0] + + xmin := (FSDWidth-f.particlebox.X)/2 + FSDParticleRadius + xmax := f.particlebox.X - FSDParticleRadius*2 + ymin := (FSDHeight-f.particlebox.Y)/2 + FSDParticleRadius + ymax := f.particlebox.Y - FSDParticleRadius*2 + + for i := 0; i < FSDParticleCount; i++ { + + p := &Particle{ + Position: gamedata.Vector{ + X: xmin + rand.Float64()*xmax, + Y: ymin + rand.Float64()*ymax, + }, + Velocity: gamedata.Vector{ + X: 0, + Y: 0, + }, + Radius: FSDParticleRadius, + } + + f.particles = append(f.particles, p) + } + +} + +func (f *FluidSimD) BoundParticle(p *Particle) { + + xmin := (FSDWidth-f.particlebox.X)/2 + p.Radius + xmax := xmin + f.particlebox.X - p.Radius*2 + + if p.Position.X > xmax { + p.Velocity.X *= -1 * FSDDamping + p.Position.X = xmax + } + + if p.Position.X < xmin { + p.Velocity.X *= -1 * FSDDamping + p.Position.X = xmin + } + + ymin := (FSDHeight-f.particlebox.Y)/2 + p.Radius + ymax := ymin + f.particlebox.Y - p.Radius*2 + + if p.Position.Y > ymax { + p.Velocity.Y *= -1 * FSDDamping + p.Position.Y = ymax + } + + if p.Position.Y < ymin { + p.Velocity.Y *= -1 * FSDDamping + p.Position.Y = ymin + } + +} + +func (f *FluidSimD) RenderQuadrants() { + clr := color.RGBA{R: 0xff, G: 0x00, B: 0x00, A: 0xff} + + quadrants := f.quadtree.GetQuadrants() + for _, quad := range quadrants { + ox := float32(quad.Position.X - quad.Dimensions.X/2) + oy := float32(quad.Position.Y - quad.Dimensions.Y/2) + vector.StrokeRect(f.Sprite, ox, oy, float32(quad.Dimensions.X), float32(quad.Dimensions.Y), 1, clr, true) + } + +} + +func (f *FluidSimD) RenderParticles() { + for _, particle := range f.particles { + x0 := particle.Position.X - particle.Radius + y0 := particle.Position.Y - particle.Radius + op := &ebiten.DrawImageOptions{} + op.GeoM.Translate(x0, y0) + f.Sprite.DrawImage(f.particlebuff, op) + } +} + +func (f *FluidSimD) RenderBox() { + x0 := (FSDWidth - f.particlebox.X) / 2 + y0 := (FSDHeight - f.particlebox.Y) / 2 + vector.StrokeRect(f.Sprite, float32(x0), float32(y0), float32(f.particlebox.X), float32(f.particlebox.Y), 2, color.White, true) +} + +func (f *FluidSimD) RebuildQuadtree() { + f.quadtree.Clear() + + for _, p := range f.particles { + if !f.quadtree.Insert(p) { + fmt.Println("quadtree insertion failed") + } + } +} + +func (f *FluidSimD) ResolveCollisionsA(particle *Particle) { + //construct search quadrant from current particle + quadrant := quadtree.Quadrant{ + Position: particle.Position, + Dimensions: particle.GetDimensions(), + } + + //find list of possible maybe collisions, we inspect those in more detail + maybes := f.quadtree.FindAll(quadrant) + + sqdist := float64(particle.Radius*particle.Radius) * 4 + + for _, p := range maybes { + if p == particle { + continue + } + + pos := p.GetPosition() + delta := gamedata.Vector{ + X: pos.X - particle.Position.X, + Y: pos.Y - particle.Position.Y, + } + dist2 := delta.X*delta.X + delta.Y*delta.Y + + if dist2 == 0 { + // Same position: pick a fallback direction to avoid NaN + delta.X = 1 + delta.Y = 0 + dist2 = 1 + } + + if dist2 < sqdist { + d := math.Sqrt(dist2) + nx, ny := delta.X/d, delta.Y/d + overlap := particle.Radius*2 - d + + pos.X += nx * overlap + pos.Y += ny * overlap + p.SetPosition(pos) + /* + newpos := gamedata.Vector{ + X: pos.X + delta.X, + Y: pos.Y + delta.Y, + } + p.SetPosition(newpos) + */ + } + } +} + +func (f *FluidSimD) ResolveCollisionsB(particle *Particle) { + //construct search quadrant from current particle + quadrant := quadtree.Quadrant{ + Position: particle.Position, + Dimensions: particle.GetDimensions(), + } + + //find list of possible maybe collisions, we inspect those in more detail + maybes := f.quadtree.FindAll(quadrant) + + sqdist := float64(particle.Radius*particle.Radius) * 4 + + for _, p := range maybes { + if p == particle { + continue + } + + pos := p.GetPosition() + delta := gamedata.Vector{ + X: pos.X - particle.Position.X, + Y: pos.Y - particle.Position.Y, + } + + dist2 := delta.X*delta.X + delta.Y*delta.Y + + if dist2 < sqdist { + d := math.Sqrt(dist2) + overlap := particle.Radius*2 - d + theta := math.Atan2(delta.Y, delta.X) + pos.X += overlap * math.Cos(theta) + pos.Y += overlap * math.Sin(theta) + p.SetPosition(pos) + + m := p.(*Particle) + m.Velocity.X *= -1 * FSDDamping + m.Velocity.Y *= -1 * FSDDamping + + } + } +} + +func (f *FluidSimD) ResolveCollisionsC(particle *Particle) { + //construct search quadrant from current particle + quadrant := quadtree.Quadrant{ + Position: particle.Position, + Dimensions: particle.GetDimensions(), + } + + //find list of possible maybe collisions, we inspect those in more detail + maybes := f.quadtree.FindAll(quadrant) + + sqdist := float64(particle.Radius*particle.Radius) * 4 + + for _, p := range maybes { + if p == particle { + continue + } + + pos := p.GetPosition() + delta := gamedata.Vector{ + X: pos.X - particle.Position.X, + Y: pos.Y - particle.Position.Y, + } + + dist2 := delta.X*delta.X + delta.Y*delta.Y + + if dist2 < sqdist { + + /* + //compute impact to this particle + deltav1 := particle.Velocity.Subtract(p.(*Particle).Velocity) + deltax1 := particle.Position.Subtract(p.(*Particle).Position) + dot1 := deltav1.DotProduct(deltax1) + mag1 := deltax1.Magnitude() * deltax1.Magnitude() + particle.Velocity = deltax1.Scale(dot1 / mag1) + + //compute impact to other particle + deltav2 := p.(*Particle).Velocity.Subtract(particle.Velocity) + deltax2 := p.(*Particle).Position.Subtract(particle.Position) + dot2 := deltav2.DotProduct(deltax2) + mag2 := deltax2.Magnitude() * deltax2.Magnitude() + p.(*Particle).Velocity = deltax2.Scale(dot2 / mag2) + */ + + dist := math.Sqrt(dist2) + s := 0.5 * (particle.Radius*2 - dist) / dist + + dnorm := delta.Scale(s) + + particle.Position = particle.Position.Subtract(dnorm) + p.(*Particle).Position = p.(*Particle).Position.Add(dnorm) + + } + } + +} + +func (f *FluidSimD) SetRenderQuads(v bool) { + f.renderquads = v +} + +func (f *FluidSimD) RenderQuads() bool { + return f.renderquads +} + +func (f *FluidSimD) SetResolveCollisions(v bool) { + f.resolvecollisions = v +} + +func (f *FluidSimD) ResolveCollisions() bool { + return f.resolvecollisions +} + +func (f *FluidSimD) NextSolver() { + f.resolveridx = (f.resolveridx + 1) % len(f.resolvers) + +} + +func (f *FluidSimD) PreviousSolver() { + f.resolveridx = f.resolveridx - 1 + if f.resolveridx < 0 { + f.resolveridx = len(f.resolvers) - 1 + } +} diff --git a/elements/mappedentity.go b/elements/mappedentity.go new file mode 100644 index 0000000..a754a34 --- /dev/null +++ b/elements/mappedentity.go @@ -0,0 +1,18 @@ +package elements + +import ( + "fluids/gamedata" + + "github.com/hajimehoshi/ebiten/v2" +) + +type MappedEntity interface { + Draw() + Update() + GetSprite() *ebiten.Image + GetDimensions() gamedata.Vector + GetPosition() gamedata.Vector + SetPosition(gamedata.Vector) + SetPaused(bool) + Paused() bool +} diff --git a/elements/mappedentitybase.go b/elements/mappedentitybase.go new file mode 100644 index 0000000..92cef87 --- /dev/null +++ b/elements/mappedentitybase.go @@ -0,0 +1,51 @@ +package elements + +import ( + "fluids/gamedata" + + "github.com/hajimehoshi/ebiten/v2" +) + +type MappedEntityBase struct { + Sprite *ebiten.Image + dimensions gamedata.Vector + position gamedata.Vector + paused bool +} + +func NewMappedEntityBase(dimensions gamedata.Vector) *MappedEntityBase { + meb := &MappedEntityBase{ + dimensions: dimensions, + } + return meb +} + +func (m *MappedEntityBase) Draw() { +} + +func (m *MappedEntityBase) Update() { +} + +func (m *MappedEntityBase) GetSprite() *ebiten.Image { + return m.Sprite +} + +func (m *MappedEntityBase) GetDimensions() gamedata.Vector { + return m.dimensions +} + +func (m *MappedEntityBase) GetPosition() gamedata.Vector { + return m.position +} + +func (m *MappedEntityBase) SetPosition(pos gamedata.Vector) { + m.position = pos +} + +func (m *MappedEntityBase) SetPaused(p bool) { + m.paused = p +} + +func (m *MappedEntityBase) Paused() bool { + return m.paused +} diff --git a/fluid/field.go b/fluid/field.go new file mode 100644 index 0000000..4026782 --- /dev/null +++ b/fluid/field.go @@ -0,0 +1,74 @@ +package fluid + +import "math" + +type CellType int + +const ( + CellTypeFluid = iota + CellTypeAir + CellTypeSolid + CellTypeMax +) + +type VelocityField struct { + Dimensions FieldVector //in meters + Numcells int + Nx, Ny int //number of cells in x, y + H float32 //field spacing, in meters + InvH float32 + U, V []float32 //field components < u(x,y), v(x,y) > + PrevU, PrevV []float32 //previous field components < u(x,y), v(x,y) > + DU, DV []float32 //weighted sums for components < u(x,y), v(x,y) > + S []float32 //fluid cell scales + CellType []float32 //type of cell +} + +func NewVelocityField(dimensions FieldVector, spacing float32) *VelocityField { + vf := &VelocityField{ + Dimensions: dimensions, + H: spacing, + InvH: 1 / spacing, + } + + vf.Initialize() + return vf +} + +func (vf *VelocityField) Initialize() { + + vf.Nx = int(math.Floor(float64(vf.Dimensions.X/vf.H))) + 1 + vf.Ny = int(math.Floor(float64(vf.Dimensions.Y/vf.H))) + 1 + + cellcount := vf.Nx * vf.Ny + + vf.Numcells = cellcount + + vf.U = make([]float32, cellcount) + vf.V = make([]float32, cellcount) + vf.PrevU = make([]float32, cellcount) + vf.PrevV = make([]float32, cellcount) + vf.DU = make([]float32, cellcount) + vf.DV = make([]float32, cellcount) + vf.S = make([]float32, cellcount) + vf.CellType = make([]float32, cellcount) + + for i := 0; i < vf.Nx; i++ { + for j := 0; j < vf.Ny; j++ { + var stype float32 = CellTypeAir //port, code corrected + //var stype float32 = 1 //10mp value, claims fluid but isn't + //var stype float32 = CellTypeSolid //port, with 10mp intent + if i == 0 || i == vf.Nx-1 || j == 0 { + stype = CellTypeFluid //port, code corrected + //stype = 0 //10mp value, claims solid but isn't + //stype = CellTypeFluid //port, with 10mp intent + } + vf.S[i*vf.Ny+j] = stype + } + } + +} + +func (vf *VelocityField) GetIndex(i, j int) int { + return i*vf.Ny + j +} diff --git a/fluid/fieldvector.go b/fluid/fieldvector.go new file mode 100644 index 0000000..a753a17 --- /dev/null +++ b/fluid/fieldvector.go @@ -0,0 +1,36 @@ +package fluid + +import "math" + +type FieldVector struct { + X float32 + Y float32 +} + +func (v FieldVector) DotProduct(p FieldVector) float32 { + return v.X*p.X + v.Y*p.Y +} + +func (v FieldVector) Magnitude() float32 { + return float32(math.Sqrt(float64(v.X*v.X + v.Y*v.Y))) +} + +func (v FieldVector) Add(p FieldVector) FieldVector { + return FieldVector{X: v.X + p.X, Y: v.Y + p.Y} +} + +func (v FieldVector) Subtract(p FieldVector) FieldVector { + return FieldVector{X: v.X - p.X, Y: v.Y - p.Y} +} + +func (v FieldVector) Scale(s float32) FieldVector { + return FieldVector{X: v.X * s, Y: v.Y * s} +} + +func (v FieldVector) Normalize() FieldVector { + mag := v.Magnitude() + return FieldVector{ + X: v.X / mag, + Y: v.Y / mag, + } +} diff --git a/fluid/fluid.go b/fluid/fluid.go new file mode 100644 index 0000000..4696f5e --- /dev/null +++ b/fluid/fluid.go @@ -0,0 +1,626 @@ +// Based on code by Matthias Müller (Ten Minute Physics), MIT License +// See licenses/TenMinutePhysics for full license text + +package fluid + +import ( + "math" +) + +const ( + FluidDefaultDensity = 1000 //kilogram per cubic meter + FluidDefaultGravity = -9.8 + FluidDefaultTimeStep = 1. / 120 + FluidDefaultParticleRadiusPercent = 0.3 + FluidDefaultAngle = 0 + FluidIncompressibilityIterations = 100 + FluidDefaultRestDensity = 0 + FluidDefaultFlipPicRatio = 0.9 + FluidDefaultOverRelaxation = 1.9 + FluidDefaultSubSteps = 2 +) + +type Fluid struct { + Particles []Particle + Field *VelocityField + particleRadius float32 + numCellParticles []int //SPATIAL HASHING: count of particles per cell + firstCellParticle []int //SPATIAL HASHING: index of first particle for given cell + cellParticleIds []int //SPATIAL HASHING: id of all particles + pNumX, pNumY int //SPATIAL HASHING: number of particle cells x,y + pNumCells int //SPATIAL HASHING: number of particle cells in total + pInvSpacing float32 //SPATIAL HASHING: one over the particle spacing + particleDensity []float32 //density of particles within each field cell + particleRestDensity float32 + fluidDensity float32 + angle float32 + dt float32 //delta time step + gravitynormal FieldVector //unit vector for gravity + stepacceleration FieldVector //step acceleration due to gravity + flipPicRatio float32 + numSubSteps int //number of simulation substeps +} + +func NewFluid(dimensions FieldVector, spacing float32) *Fluid { + f := &Fluid{ + Field: NewVelocityField(dimensions, spacing), + dt: FluidDefaultTimeStep, + fluidDensity: FluidDefaultDensity, + particleRestDensity: FluidDefaultRestDensity, + flipPicRatio: FluidDefaultFlipPicRatio, + numSubSteps: FluidDefaultSubSteps, + } + + f.Initialize() + return f +} + +func (f *Fluid) Initialize() { + f.InitializeParticles() + + f.SetAngle(FluidDefaultAngle) + f.particleDensity = make([]float32, f.Field.Nx*f.Field.Ny) +} + +func (f *Fluid) InitializeParticles() { + //compute particle radius based on our field spacing and radius scaling + f.particleRadius = FluidDefaultParticleRadiusPercent * f.Field.H + + //we want to prepare a hexagonally packed grid of particles within our fluid field + dx := float32(2.0 * f.particleRadius) + dy := float32(math.Sqrt(3.) / 2. * float64(dx)) + + //number of particles that fit within our specified dimensions + numPx := int(math.Floor(float64(f.Field.Dimensions.X-2*f.Field.H-2*f.particleRadius) / float64(dx))) + numPy := int(math.Floor(float64(f.Field.Dimensions.Y/2-2*f.Field.H-2*f.particleRadius) / float64(dy))) + + particlecount := numPx * numPy + f.Particles = make([]Particle, particlecount) + + //now prepare our hashing structures, purely for collision adjustment + //first let's compute some boundaries + f.pInvSpacing = 1 / (2.2 * f.particleRadius) + f.pNumX = int(math.Floor(float64(f.Field.Dimensions.X*f.pInvSpacing))) + 1 + f.pNumY = int(math.Floor(float64(f.Field.Dimensions.Y*f.pInvSpacing))) + 1 + f.pNumCells = f.pNumX * f.pNumY + + //now use boundaries to construct our spatial hashing containers + f.numCellParticles = make([]int, f.pNumCells) + f.firstCellParticle = make([]int, f.pNumCells+1) + f.cellParticleIds = make([]int, int(particlecount)) + + //create particles + var p int = 0 + for i := 0; i < numPx; i++ { + for j := 0; j < numPy; j++ { + + var rowAdjust float32 = 0 + if j%2 == 1 { + rowAdjust = f.particleRadius + } + + f.Particles[p].Position.X = f.Field.H + f.particleRadius + dx*float32(i) + rowAdjust + f.Particles[p].Position.Y = f.Field.H + f.particleRadius + dy*float32(j) + + p++ + } + } +} + +func (f *Fluid) Step() { + + for range f.numSubSteps { + f.IntegrateParticles() + f.SeparateParticles() + f.HandleParticleCollisions() + f.TransferVelocities(true) + f.UpdateParticleDensity() + f.SolveIncompressibility() + f.TransferVelocities(false) + } + //f.UpdateParticleColors() +} + +func (f *Fluid) RecalculateStepAcceleration() { + scale := FluidDefaultGravity * f.dt + f.stepacceleration = FieldVector{ + X: scale * f.gravitynormal.X, + Y: scale * f.gravitynormal.Y, + } +} + +// change the time step +func (f *Fluid) SetTimeStep(dt float32) { + f.dt = dt + f.RecalculateStepAcceleration() +} + +// set new fluid angle, recompute gravity normal, initiate step accel recomputation +func (f *Fluid) SetAngle(angle float32) { + f.angle = angle + sin, cos := math.Sincos(float64(angle)) + f.gravitynormal = FieldVector{ + X: -float32(sin), + Y: float32(cos), + } + f.RecalculateStepAcceleration() +} + +func (f *Fluid) IntegrateParticles() { + ax, ay := f.stepacceleration.X, f.stepacceleration.Y + dt := f.dt + for i := range f.Particles { + p := &f.Particles[i] + p.Velocity.X += ax + p.Velocity.Y += ay + p.Position.X += p.Velocity.X * dt + p.Position.Y += p.Velocity.Y * dt + } +} + +// returns adjusted value between start and end, whichever is closer +// unless value is between, in which case it returns value +func clamp(value, start, end int) int { + if value < start { + return start + } + + if value < end { + return value + } + + return end +} + +func clamp32(value, start, end float32) float32 { + if value < start { + return start + } + + if value < end { + return value + } + + return end +} + +func (f *Fluid) PerformParticleSpatialHash() { + + clear(f.numCellParticles) + + //find spatcial hash cell for this particle's position, increment count there + for i := range f.Particles { + p := &f.Particles[i] + xi := clamp(int(p.Position.X*f.pInvSpacing), 0, f.pNumX-1) + yi := clamp(int(p.Position.Y*f.pInvSpacing), 0, f.pNumY-1) + cellnumber := xi*f.pNumY + yi + f.numCellParticles[cellnumber]++ + } + + //build up firstCellParticle bucket + var sum int = 0 + for i := range f.pNumCells { + sum += f.numCellParticles[i] + f.firstCellParticle[i] = sum + } + f.firstCellParticle[f.pNumCells] = sum + + //next work backwards on it to assign particle ids, adjusting the + //firstCellParticle bucket as we go + for i := range f.Particles { + p := &f.Particles[i] + xi := clamp(int(p.Position.X*f.pInvSpacing), 0, f.pNumX-1) + yi := clamp(int(p.Position.Y*f.pInvSpacing), 0, f.pNumY-1) + cellnumber := xi*f.pNumY + yi + f.firstCellParticle[cellnumber]-- + f.cellParticleIds[f.firstCellParticle[cellnumber]] = i + } +} + +func (f *Fluid) SeparateParticles() { + + //setup + minDist := 2 * f.particleRadius + minDist2 := minDist * minDist + + //build our spatial hashing table + f.PerformParticleSpatialHash() + + //perform broad phase search + for i := range f.Particles { + + p := &f.Particles[i] + + pxi := int(p.Position.X * f.pInvSpacing) + pyi := int(p.Position.Y * f.pInvSpacing) + x0 := max(pxi-1, 0) + y0 := max(pyi-1, 0) + x1 := min(pxi+1, f.pNumX) + y1 := min(pyi+1, f.pNumY) + + for xi := x0; xi < x1; xi++ { + for yi := y0; yi < y1; yi++ { + //find current cell + cell := xi*f.pNumY + yi + + //set start and ending particles to check between this and next cell + first := f.firstCellParticle[cell] + last := f.firstCellParticle[cell+1] + + //iterate over these candidates + for j := first; j < last; j++ { + //skip self + id := f.cellParticleIds[j] + if id == i { + continue + } + + //compute delta + dx := f.Particles[id].Position.X - p.Position.X + dy := f.Particles[id].Position.Y - p.Position.Y + dist2 := dx*dx + dy*dy + + //rule out non-collisions, ignore overlaps (we screwed) + if (dist2 == 0) || (dist2 > minDist2) { + continue + } + + //need to use more expensive sqrt now to find our penetration + //then compute a scalar for vector normalization, applied + //equally in half to both particles + dist := float32(math.Sqrt(float64(dist2))) + penetration := (minDist - dist) + scale := penetration / dist + dx *= scale / 2 + dy *= scale / 2 + + p.Position.X -= dx + p.Position.Y -= dy + f.Particles[id].Position.X += dx + f.Particles[id].Position.Y += dy + } + } + } + } +} + +func (f *Fluid) HandleParticleCollisions() { + //setup + /* + minDist := 2 * f.particleRadius + minDist2 := minDist * minDist + */ + + minX := f.particleRadius + maxX := float32(f.Field.Nx-1)*f.Field.H - f.particleRadius + minY := f.particleRadius + maxY := float32(f.Field.Ny-1)*f.Field.H - f.particleRadius + + for i := range f.Particles { + p := &f.Particles[i] + + if p.Position.X < minX { + p.Position.X = minX + p.Velocity.X = 0 + } + + if p.Position.X > maxX { + p.Position.X = maxX + p.Velocity.X = 0 + } + + if p.Position.Y < minY { + p.Position.Y = minY + p.Velocity.Y = 0 + } + + if p.Position.Y > maxY { + p.Position.Y = maxY + p.Velocity.Y = 0 + } + + } +} + +func (f *Fluid) TransferVelocities(togrid bool) { + + //setup + n := f.Field.Ny + h1 := f.Field.InvH + h2 := 0.5 * f.Field.H + + if togrid { + copy(f.Field.PrevU, f.Field.U) + copy(f.Field.PrevV, f.Field.V) + clear(f.Field.DU) + clear(f.Field.DV) + clear(f.Field.U) + clear(f.Field.V) + + for i := 0; i < f.Field.Numcells; i++ { + if f.Field.S[i] == 0 { + f.Field.CellType[i] = CellTypeSolid + } else { + f.Field.CellType[i] = CellTypeAir + } + } + + for i := range f.Particles { + x := f.Particles[i].Position.X + y := f.Particles[i].Position.Y + xi := clamp(int(math.Floor(float64(x*h1))), 0, f.Field.Nx-1) + yi := clamp(int(math.Floor(float64(y*h1))), 0, f.Field.Ny-1) + cellnumber := xi*n + yi + if f.Field.CellType[cellnumber] == CellTypeAir { + f.Field.CellType[cellnumber] = CellTypeFluid + } + } + } + + for component := 0; component < 2; component++ { + var dx, dy float32 + var fl, prevF, d *[]float32 + + if component == 0 { + dx = 0 + dy = h2 + fl = &f.Field.U + prevF = &f.Field.PrevU + d = &f.Field.DU + } else { + dx = h2 + dy = 0 + fl = &f.Field.V + prevF = &f.Field.PrevV + d = &f.Field.DV + } + + for i := range f.Particles { + p := &f.Particles[i] + x := p.Position.X + y := p.Position.Y + + x = clamp32(x, f.Field.H, float32((f.Field.Nx-1))*f.Field.H) + y = clamp32(y, f.Field.H, float32((f.Field.Ny-1))*f.Field.H) + + //find cell components neighbouring x,y + x0 := min(int(math.Floor(float64((x-dx)*h1))), f.Field.Nx-2) + x1 := min(x0+1, f.Field.Nx-2) + y0 := min(int(math.Floor(float64((y-dy)*h1))), f.Field.Ny-2) + y1 := min(y0+1, f.Field.Ny-2) + + //compute scales + tx := ((x - dx) - float32(x0)*f.Field.H) * h1 + ty := ((y - dy) - float32(y0)*f.Field.H) * h1 + sx := 1.0 - tx + sy := 1.0 - ty + + //compute weights + d0 := sx * sy + d1 := tx * sy + d2 := tx * ty + d3 := sx * ty + + //cell indexes + nr0 := x0*n + y0 + nr1 := x1*n + y0 + nr2 := x1*n + y1 + nr3 := x0*n + y1 + + if togrid { + //extract velocity for this component + var pv float32 // = this.particleVel[2 * i + component]; + if component == 0 { + pv = p.Velocity.X + } else { + pv = p.Velocity.Y + } + + //to our velocity component we will add the contribution of + //the particle velocity scaled according to its weight + (*fl)[nr0] += pv * d0 + (*fl)[nr1] += pv * d1 + (*fl)[nr2] += pv * d2 + (*fl)[nr3] += pv * d3 + + //we also keep a running total of all the weights (d is r in the notes) + (*d)[nr0] += d0 + (*d)[nr1] += d1 + (*d)[nr2] += d2 + (*d)[nr3] += d3 + } else { + + var offset int + var v float32 + if component == 0 { + offset = n + v = p.Velocity.X + } else { + offset = 1 + v = p.Velocity.Y + } + + var valid0, valid1, valid2, valid3 float32 = 0., 0., 0., 0. + + if f.Field.CellType[nr0] != CellTypeAir || f.Field.CellType[nr0-offset] != CellTypeAir { + valid0 = 1 + } + + if f.Field.CellType[nr1] != CellTypeAir || f.Field.CellType[nr1-offset] != CellTypeAir { + valid1 = 1 + } + + if f.Field.CellType[nr2] != CellTypeAir || f.Field.CellType[nr2-offset] != CellTypeAir { + valid2 = 1 + } + + if f.Field.CellType[nr3] != CellTypeAir || f.Field.CellType[nr3-offset] != CellTypeAir { + valid3 = 1 + } + + //weights + wsum := valid0*d0 + valid1*d1 + valid2*d2 + valid3*d3 + + if wsum > 0.0 { + picV := (valid0*d0*(*fl)[nr0] + valid1*d1*(*fl)[nr1] + valid2*d2*(*fl)[nr2] + valid3*d3*(*fl)[nr3]) / wsum + corr := (valid0*d0*((*fl)[nr0]-(*prevF)[nr0]) + valid1*d1*((*fl)[nr1]-(*prevF)[nr1]) + + valid2*d2*((*fl)[nr2]-(*prevF)[nr2]) + valid3*d3*((*fl)[nr3]-(*prevF)[nr3])) / wsum + flipV := v + corr + vNew := (1-f.flipPicRatio)*picV + f.flipPicRatio*flipV + + if component == 0 { + p.Velocity.X = vNew + } else { + p.Velocity.Y = vNew + } + } + } + } + + if togrid { + //finally, for all q's, we divide by the weighted sum of that cell + for i := 0; i < len(*fl); i++ { + if (*d)[i] > 0 { + (*fl)[i] = (*fl)[i] / (*d)[i] + } + } + + //next we restore cells are solid to their previous values + for i := 0; i < f.Field.Nx; i++ { + for j := 0; j < f.Field.Ny; j++ { + idx := i*n + j + issolid := f.Field.CellType[idx] == CellTypeSolid + + if issolid || (i > 0 && f.Field.CellType[(i-1)*n+j] == CellTypeSolid) { + f.Field.U[idx] = f.Field.PrevU[idx] + } + + if issolid || (j > 0 && f.Field.CellType[i*n+j-1] == CellTypeSolid) { + f.Field.V[idx] = f.Field.PrevV[idx] + } + + } + } + } + + } + +} + +func (f *Fluid) UpdateParticleDensity() { + n := f.Field.Ny + h := f.Field.H + h1 := f.Field.InvH + h2 := 0.5 * h + + clear(f.particleDensity) + + for i := range f.Particles { + p := &f.Particles[i] + x := p.Position.X + y := p.Position.Y + + x = clamp32(x, h, float32(f.Field.Nx-1)*h) + y = clamp32(y, h, float32(f.Field.Ny-1)*h) + + x0 := int(math.Floor(float64((x - h2) * h1))) + x1 := min(x0+1, f.Field.Nx-2) + + y0 := int(math.Floor(float64((y - h2) * h1))) + y1 := min(y0+1, f.Field.Ny-2) + + tx := ((x - h2) - float32(x0)*h) * h1 + ty := ((y - h2) - float32(y0)*h) * h1 + + sx := 1.0 - tx + sy := 1.0 - ty + + if x0 < f.Field.Nx && y0 < f.Field.Ny { + f.particleDensity[x0*n+y0] += sx * sy + } + if x1 < f.Field.Nx && y0 < f.Field.Ny { + f.particleDensity[x1*n+y0] += tx * sy + } + if x1 < f.Field.Nx && y1 < f.Field.Ny { + f.particleDensity[x1*n+y1] += tx * ty + } + if x0 < f.Field.Nx && y1 < f.Field.Ny { + f.particleDensity[x0*n+y1] += sx * ty + } + } + + if f.particleRestDensity == 0.0 { + var sum float32 = 0.0 + numFluidCells := 0 + + for i := 0; i < f.Field.Nx*f.Field.Ny; i++ { + if f.Field.CellType[i] == CellTypeFluid { + sum += f.particleDensity[i] + numFluidCells++ + } + } + + if numFluidCells > 0 { + f.particleRestDensity = sum / float32(numFluidCells) + } + } +} + +func (f *Fluid) SolveIncompressibility() { + + //setup + copy(f.Field.PrevU, f.Field.U) + copy(f.Field.PrevV, f.Field.V) + + n := f.Field.Ny + //cp := f.fluidDensity * f.Field.H / FluidDefaultTimeStep + + for iter := 0; iter < FluidIncompressibilityIterations; iter++ { + for i := 1; i < f.Field.Nx-1; i++ { + for j := 1; j < f.Field.Ny-1; j++ { + + if f.Field.CellType[i*n+j] != CellTypeFluid { + continue + } + + //compute indexes of surrounding cells + center := i*n + j + left := (i-1)*n + j + right := (i+1)*n + j + bottom := i*n + j - 1 + top := i*n + j + 1 + + sx0 := f.Field.S[left] + sx1 := f.Field.S[right] + sy0 := f.Field.S[bottom] + sy1 := f.Field.S[top] + s := sx0 + sx1 + sy0 + sy1 + if s == 0 { + continue + } + + //divergence -> we want this to be zero + div := f.Field.U[right] - f.Field.U[center] + f.Field.V[top] - f.Field.V[center] + + if f.particleRestDensity > 0.0 { + var k float32 = 1.0 + var compression float32 = f.particleDensity[i*n+j] - f.particleRestDensity + if compression > 0.0 { + div = div - k*compression + } + } + + p := -div / s + p *= FluidDefaultOverRelaxation //overRelaxation + //f.Field.Pressure[center] += cp * p + + f.Field.U[center] -= sx0 * p + f.Field.U[right] += sx1 * p + f.Field.V[center] -= sy0 * p + f.Field.V[top] += sy1 * p + + } + } + } + +} diff --git a/fluid/particle.go b/fluid/particle.go new file mode 100644 index 0000000..5156a02 --- /dev/null +++ b/fluid/particle.go @@ -0,0 +1,6 @@ +package fluid + +type Particle struct { + Position FieldVector + Velocity FieldVector +} diff --git a/game/game.go b/game/game.go index 658ce9d..de36e66 100644 --- a/game/game.go +++ b/game/game.go @@ -3,93 +3,72 @@ package game import ( "fluids/elements" "fluids/gamedata" - "fluids/quadtree" "fmt" - "image/color" "math" - "math/rand" "github.com/hajimehoshi/ebiten/v2" + "github.com/hajimehoshi/ebiten/v2/ebitenutil" "github.com/hajimehoshi/ebiten/v2/inpututil" - "github.com/hajimehoshi/ebiten/v2/vector" ) const ( - GameWidth = 640 - GameHeight = 360 - GameParticleCount = 2000 - GameGravity = 2 - GameParticleRadius = 2.5 - GameDamping = .7 - GameDeltaTimeStep = 0.5 - GameInfluenceRadius = 30 + GameWidth = 640 + GameHeight = 360 + GameFSDW = 200 + GameFSDH = 100 ) type Game struct { - particles []*elements.Particle - cycle int - particlebox *gamedata.Vector - particlebuff *ebiten.Image - quadtree *quadtree.Quadtree - collisionquad quadtree.Quadrant - paused bool - renderquads bool - resolvecollisions bool - resolvers []func(particle *elements.Particle) - resolveridx int - alertbox *elements.Alert + //control data + cycle int + paused bool + fluidsimidx int + + //key elements + fluidsimd *elements.FluidSimD + fluidsim10 []*elements.FluidSim10 + alertbox *elements.Alert + + //cache elements + fluidsim10width float64 + fluidsim10height float64 + + //other + fluidsim10angle float64 + mdown bool + mdx, mdy int } func NewGame() *Game { g := &Game{ - particlebuff: ebiten.NewImage(GameParticleRadius*2, GameParticleRadius*2), - paused: false, - renderquads: false, - resolvecollisions: false, - resolveridx: 0, - alertbox: elements.NewAlert(), + paused: false, + fluidsimidx: 0, + alertbox: elements.NewAlert(), + fluidsimd: elements.NewFluidSimD(), + //fluidsim10: elements.NewFluidSim10(gamedata.Vector{X: GameFSDW, Y: GameFSDH}), + //fluidsimgpt: elements.NewFlipFluidEntity(640, 480, 2, 1, 100), + + //fluidsim10: elements.NewFluidSim10(), } - g.particlebox = &gamedata.Vector{ - X: GameWidth - 50, - Y: GameHeight - 50, - } - - quad := quadtree.Quadrant{ - Position: gamedata.Vector{X: GameWidth / 2, Y: GameHeight / 2}, - Dimensions: gamedata.Vector{X: GameWidth, Y: GameHeight}, - } - g.quadtree = quadtree.New(quad, 0) - - vector.FillCircle(g.particlebuff, GameParticleRadius, GameParticleRadius, GameParticleRadius, color.White, true) - - g.collisionquad = quadtree.Quadrant{ - Position: gamedata.Vector{ - X: 0, - Y: 0, - }, - Dimensions: gamedata.Vector{ - X: GameParticleRadius * 2, - Y: GameParticleRadius * 2, - }, - } - - g.resolvers = append(g.resolvers, g.ResolveCollisionsA) - g.resolvers = append(g.resolvers, g.ResolveCollisionsB) - - //g.InitializeColliders() - g.InitializeParticles() - + g.Initialize() return g } func (g *Game) Update() error { g.ParseInputs() - g.RebuildQuadtree() if !g.paused { - g.UpdateParticles() + switch g.fluidsimidx { + case 0: + g.fluidsimd.Update() + case 1: + //g.fluidsimgpt.Update() + g.UpdateFluidsim10() + default: + break + } } g.cycle++ @@ -97,13 +76,22 @@ func (g *Game) Update() error { return nil } +func (g *Game) UpdateFluidsim10() { + for _, sim := range g.fluidsim10 { + sim.Update() + } +} + func (g *Game) Draw(screen *ebiten.Image) { screen.Clear() - g.RenderParticles(screen) - g.RenderBox(screen) - if g.renderquads { - g.RenderQuadrants(screen) + switch g.fluidsimidx { + case 0: + g.RenderFluidSimD(screen) + case 1: + g.RenderFluidSim10(screen) + default: + break } if g.paused { @@ -113,314 +101,161 @@ func (g *Game) Draw(screen *ebiten.Image) { } } +func (g *Game) RenderFluidSimD(img *ebiten.Image) { + g.fluidsimd.Draw() + pos := g.fluidsimd.GetPosition() + op := &ebiten.DrawImageOptions{} + op.GeoM.Translate(pos.X, pos.Y) + img.DrawImage(g.fluidsimd.GetSprite(), op) +} + +func (g *Game) RenderFluidSim10(img *ebiten.Image) { + + for _, sim := range g.fluidsim10 { + sim.Draw() + + pos := sim.GetPosition() + op := &ebiten.DrawImageOptions{} + op.GeoM.Translate(-g.fluidsim10width/2, -g.fluidsim10height/2) + op.GeoM.Scale(1, -1) + op.GeoM.Rotate(g.fluidsim10angle) + // op.GeoM.Translate(g.fluidsim10width/2, g.fluidsim10height/2) + op.GeoM.Translate(pos.X, pos.Y) + img.DrawImage(sim.GetSprite(), op) + } + + //debug info + x, y := ebiten.CursorPosition() + deg := g.fluidsim10angle * 180 / math.Pi + + str := fmt.Sprintf("Mouse (x: %d, y: %d) Origin (x: %d, y: %d) Angle (%f)", x, y, GameWidth/2, GameHeight/2, deg) + ebitenutil.DebugPrint(img, str) +} + func (g *Game) Layout(x, y int) (int, int) { return GameWidth, GameHeight } -func (g *Game) RenderParticles(img *ebiten.Image) { - // mx, my := ebiten.CursorPosition() - //clr := color.RGBA{R: 0xff, G: 0x00, B: 0x00, A: 0xff} - - for _, particle := range g.particles { - - x0 := particle.Position.X - particle.Radius - y0 := particle.Position.Y - particle.Radius - - //redness := float32(particle.Position.Y / g.particlebox.Y) - //blueness := 1 - redness - - op := &ebiten.DrawImageOptions{} - op.GeoM.Translate(x0, y0) - //op.ColorScale.Scale(redness, 0, blueness, 1) - img.DrawImage(g.particlebuff, op) - - //vector.FillCircle(img, float32(particle.Position.X), float32(particle.Position.Y), float32(particle.Radius), color.White, true) - // vector.StrokeCircle(img, float32(particle.Position.X), float32(particle.Position.Y), GameInfluenceRadius, 1, clr, true) - // vector.StrokeLine(img, float32(mx), float32(my), float32(particle.Position.X), float32(particle.Position.Y), 2, clr, true) - } -} - -func (g *Game) RenderBox(img *ebiten.Image) { - x0 := (GameWidth - g.particlebox.X) / 2 - y0 := (GameHeight - g.particlebox.Y) / 2 - vector.StrokeRect(img, float32(x0), float32(y0), float32(g.particlebox.X), float32(g.particlebox.Y), 2, color.White, true) -} - -func (g *Game) InitializeColliders() { - /* - //box container - box := &colliders.BaseRectCollider{} - box.SetDimensions(gamedata.Vector{ - X: GameWidth - 50, - Y: GameHeight - 50, - }) - box.SetPosition(gamedata.Vector{ - X: GameWidth / 2, - Y: GameHeight / 2, - }) - box.SetContainer(true) - g.rectColliders = append(g.rectColliders, box) - */ -} - -func (g *Game) InitializeParticles() { - - g.particles = g.particles[:0] - - xmin := (GameWidth-g.particlebox.X)/2 + GameParticleRadius - xmax := g.particlebox.X - GameParticleRadius*2 - ymin := (GameHeight-g.particlebox.Y)/2 + GameParticleRadius - ymax := g.particlebox.Y - GameParticleRadius*2 - - for i := 0; i < GameParticleCount; i++ { - - p := &elements.Particle{ - Position: gamedata.Vector{ - X: xmin + rand.Float64()*xmax, - Y: ymin + rand.Float64()*ymax, - }, - Velocity: gamedata.Vector{ - X: 0, - Y: 0, - }, - Radius: GameParticleRadius, - } - - g.particles = append(g.particles, p) - } - -} - -func (g *Game) UpdateParticles() { - - dt := GameDeltaTimeStep - - mx, my := ebiten.CursorPosition() - mpos := gamedata.Vector{X: float64(mx), Y: float64(my)} - maxdeflect := 40 * GameDeltaTimeStep * GameGravity - - for _, particle := range g.particles { - particle.Velocity.Y += GameGravity * dt - - //if inpututil.IsMouseButtonJustPressed(ebiten.MouseButtonLeft) { - - if ebiten.IsMouseButtonPressed(ebiten.MouseButtonLeft) { - delta := gamedata.Vector{ - X: mpos.X - particle.Position.X, - Y: mpos.Y - particle.Position.Y, - } - - dist := math.Sqrt(delta.X*delta.X + delta.Y*delta.Y) - theta := math.Atan2(delta.Y, delta.X) - - if dist < GameInfluenceRadius { - - dx := dist * math.Cos(theta) - dy := dist * math.Sin(theta) - - if dx != 0 { - gainx := (-1./GameInfluenceRadius)*math.Abs(dx) + 1. - particle.Velocity.X += 20 * gainx * -1 * math.Copysign(1, dx) - } - - if dy != 0 { - gainy := (-1./GameInfluenceRadius)*math.Abs(dy) + 1. - particle.Velocity.Y += maxdeflect * gainy * -1 * math.Copysign(1, dy) - } - } - - } - - particle.Position.X += particle.Velocity.X * dt - particle.Position.Y += particle.Velocity.Y * dt - - if g.resolvecollisions { - g.resolvers[g.resolveridx](particle) - } - } - - for _, p := range g.particles { - g.BoundParticle(p) - } - -} - -func (g *Game) BoundParticle(p *elements.Particle) { - - xmin := (GameWidth-g.particlebox.X)/2 + p.Radius - xmax := xmin + g.particlebox.X - p.Radius*2 - - if p.Position.X > xmax { - p.Velocity.X *= -1 * GameDamping - p.Position.X = xmax - } - - if p.Position.X < xmin { - p.Velocity.X *= -1 * GameDamping - p.Position.X = xmin - } - - ymin := (GameHeight-g.particlebox.Y)/2 + p.Radius - ymax := ymin + g.particlebox.Y - p.Radius*2 - - if p.Position.Y > ymax { - p.Velocity.Y *= -1 * GameDamping - p.Position.Y = ymax - } - - if p.Position.Y < ymin { - p.Velocity.Y *= -1 * GameDamping - p.Position.Y = ymin - } - -} - -func (g *Game) RenderQuadrants(img *ebiten.Image) { - clr := color.RGBA{R: 0xff, G: 0x00, B: 0x00, A: 0xff} - - quadrants := g.quadtree.GetQuadrants() - for _, quad := range quadrants { - ox := float32(quad.Position.X - quad.Dimensions.X/2) - oy := float32(quad.Position.Y - quad.Dimensions.Y/2) - vector.StrokeRect(img, ox, oy, float32(quad.Dimensions.X), float32(quad.Dimensions.Y), 1, clr, true) - } - -} - func (g *Game) ParseInputs() { - //refresh particles - if inpututil.IsKeyJustPressed(ebiten.KeyR) { - g.InitializeParticles() + //simulation specific updates + switch g.fluidsimidx { + case 0: + g.ManageFluidSimDInputs() + case 1: + //g.ManageFluidSimGPTInputs() + g.ManageFluidSim10Inputs() + default: + break } - //pause simulation + //common updates if inpututil.IsKeyJustPressed(ebiten.KeyP) { g.paused = !g.paused g.alertbox.SetText("PAUSED") g.alertbox.Draw() } + //swap fluid simulations + if inpututil.IsKeyJustPressed(ebiten.KeyPageUp) { + g.fluidsimidx = (g.fluidsimidx + 1) % 2 + } + + if inpututil.IsKeyJustPressed(ebiten.KeyPageDown) { + g.fluidsimidx = g.fluidsimidx - 1 + if g.fluidsimidx < 0 { + g.fluidsimidx = 1 + } + } + +} + +func (g *Game) ManageFluidSimDInputs() { + //refresh particles + if inpututil.IsKeyJustPressed(ebiten.KeyR) { + g.fluidsimd.InitializeParticles() + } + + //pause simulation + if inpututil.IsKeyJustPressed(ebiten.KeyP) { + g.fluidsimd.SetPaused(!g.fluidsimd.Paused()) + } + //show quadtree quadrants if inpututil.IsKeyJustPressed(ebiten.KeyQ) { - g.renderquads = !g.renderquads + g.fluidsimd.SetRenderQuads(!g.fluidsimd.RenderQuads()) } //enable collision resolution if inpututil.IsKeyJustPressed(ebiten.KeyC) { - g.resolvecollisions = !g.resolvecollisions + g.fluidsimd.SetResolveCollisions(!g.fluidsimd.ResolveCollisions()) } //switch between collision resolvers if inpututil.IsKeyJustPressed(ebiten.KeyLeft) { - g.resolveridx = g.resolveridx - 1 - if g.resolveridx < 0 { - g.resolveridx = len(g.resolvers) - 1 - } + g.fluidsimd.PreviousSolver() } if inpututil.IsKeyJustPressed(ebiten.KeyRight) { - g.resolveridx = (g.resolveridx + 1) % len(g.resolvers) + g.fluidsimd.NextSolver() } } -func (g *Game) RebuildQuadtree() { - g.quadtree.Clear() - - for _, p := range g.particles { - if !g.quadtree.Insert(p) { - fmt.Println("quadtree insertion failed") +func (g *Game) ManageFluidSim10Inputs() { + //refresh particles + if inpututil.IsKeyJustPressed(ebiten.KeyR) { + for _, sim := range g.fluidsim10 { + sim.Initialize() + g.fluidsim10angle = 0 } } + + if inpututil.IsMouseButtonJustPressed(ebiten.MouseButtonLeft) { + g.mdown = true + g.mdx, g.mdy = ebiten.CursorPosition() + } + + if ebiten.IsMouseButtonPressed(ebiten.MouseButtonLeft) { + if g.mdown { + mx, my := ebiten.CursorPosition() + dx := float64(mx) - GameWidth/2 //g.fluidsim10.GetPosition().X + dy := float64(my) - GameHeight/2 //g.fluidsim10.GetPosition().Y + angle := math.Atan2(dy, dx) + g.fluidsim10angle = angle + + for _, sim := range g.fluidsim10 { + sim.SetAngle(angle) + } + } + } + + if inpututil.IsMouseButtonJustReleased(ebiten.MouseButtonLeft) { + g.mdown = false + } + } -func (g *Game) ResolveCollisionsA(particle *elements.Particle) { - //construct search quadrant from current particle - quadrant := quadtree.Quadrant{ - Position: particle.Position, - Dimensions: particle.GetDimensions(), - } +func (g *Game) ManageFluidSimGPTInputs() { - //find list of possible maybe collisions, we inspect those in more detail - maybes := g.quadtree.FindAll(quadrant) - - sqdist := float64(particle.Radius*particle.Radius) * 4 - - for _, p := range maybes { - if p == particle { - continue - } - - pos := p.GetPosition() - delta := gamedata.Vector{ - X: pos.X - particle.Position.X, - Y: pos.Y - particle.Position.Y, - } - dist2 := delta.X*delta.X + delta.Y*delta.Y - - if dist2 == 0 { - // Same position: pick a fallback direction to avoid NaN - delta.X = 1 - delta.Y = 0 - dist2 = 1 - } - - if dist2 < sqdist { - d := math.Sqrt(dist2) - nx, ny := delta.X/d, delta.Y/d - overlap := particle.Radius*2 - d - - pos.X += nx * overlap - pos.Y += ny * overlap - p.SetPosition(pos) - /* - newpos := gamedata.Vector{ - X: pos.X + delta.X, - Y: pos.Y + delta.Y, - } - p.SetPosition(newpos) - */ - } - } } -func (g *Game) ResolveCollisionsB(particle *elements.Particle) { - //construct search quadrant from current particle - quadrant := quadtree.Quadrant{ - Position: particle.Position, - Dimensions: particle.GetDimensions(), +func (g *Game) Initialize() { + + g.fluidsim10 = append(g.fluidsim10, elements.NewFluidSim10()) + g.fluidsim10 = append(g.fluidsim10, elements.NewFluidSim10()) + + g.fluidsim10width = float64(g.fluidsim10[0].GetSprite().Bounds().Dx()) + g.fluidsim10height = float64(g.fluidsim10[0].GetSprite().Bounds().Dy()) + + x0 := float64(GameWidth / (len(g.fluidsim10) + 1)) + + for i, sim := range g.fluidsim10 { + + pos := gamedata.Vector{ + X: x0 * float64(i+1), + Y: GameHeight / 2., + } + sim.SetPosition(pos) } - //find list of possible maybe collisions, we inspect those in more detail - maybes := g.quadtree.FindAll(quadrant) - - sqdist := float64(particle.Radius*particle.Radius) * 4 - - for _, p := range maybes { - if p == particle { - continue - } - - pos := p.GetPosition() - delta := gamedata.Vector{ - X: pos.X - particle.Position.X, - Y: pos.Y - particle.Position.Y, - } - - dist2 := delta.X*delta.X + delta.Y*delta.Y - - if dist2 < sqdist { - d := math.Sqrt(dist2) - overlap := particle.Radius*2 - d - theta := math.Atan2(delta.Y, delta.X) - pos.X += overlap * math.Cos(theta) - pos.Y += overlap * math.Sin(theta) - p.SetPosition(pos) - - m := p.(*elements.Particle) - m.Velocity.X *= -1 * GameDamping - m.Velocity.Y *= -1 * GameDamping - - } - } } diff --git a/gamedata/vector.go b/gamedata/vector.go index becdff4..fcb986f 100644 --- a/gamedata/vector.go +++ b/gamedata/vector.go @@ -1,6 +1,28 @@ package gamedata +import "math" + type Vector struct { X float64 Y float64 } + +func (v Vector) DotProduct(p Vector) float64 { + return v.X*p.X + v.Y + p.Y +} + +func (v Vector) Magnitude() float64 { + return math.Sqrt(v.X*v.X + v.Y*v.Y) +} + +func (v Vector) Add(p Vector) Vector { + return Vector{X: v.X + p.X, Y: v.Y + p.Y} +} + +func (v Vector) Subtract(p Vector) Vector { + return Vector{X: v.X - p.X, Y: v.Y - p.Y} +} + +func (v Vector) Scale(s float64) Vector { + return Vector{X: v.X * s, Y: v.Y * s} +} diff --git a/licenses/TenMinutePhysics b/licenses/TenMinutePhysics new file mode 100644 index 0000000..fd1604d --- /dev/null +++ b/licenses/TenMinutePhysics @@ -0,0 +1,11 @@ +Copyright 2022 Matthias Müller - Ten Minute Physics, +www.youtube.com/c/TenMinutePhysics +www.matthiasMueller.info/tenMinutePhysics + +MIT License + +Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. \ No newline at end of file