fluids/elements/fluidsimd.go

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
	}
}