Barnes-Hut

From a question I posted, a gentleman suggested that I ask on the forums.

http://answers.unity3d.com/questions/882719/barnes-hut.html

Ok I found http://beltoforion.de/barnes-hut/barnes-hut_en.html which actually is quite helpful. I am reading up on quad trees and im just wondering if unity has such a data structure?

 

Ok I honestly at this point have no idea how to implement this. I understand how it works however coding it seems to be a major problem (the quadtree). The force of gravity and center of mass part is very simple which I already done but brute forcing is causing it to be pretty inefficient (O(N^2))

 

How many objects are you aiming for?

 

Well if I have 1000 stars that is at most about 10000 planets and 50000 moons therefore, 61000 objects. This number would be much smaller of course but Barnes-hut should be able to easily handle that amount.

I think the confusion comes from the quad tree itself but I can't find any code or good articles that I can study. They all seem to be conceptual or I'm not a good google searcher... I don't want any exact implementation just some detailed articles or tutorials so I can learn the material by filling in the blank (for example pseudo code)

 

http://www.flipcode.com/archives/Octree_Implementation.shtml
You can easily convert this to c#

Does it have to work in real time? I suspect even with Barnes-hut it would be possible in real time for 61000 particles.

 

Ok I looked at a java implementation and rewrote the BarnesHut to work with my main NBody gravity manager for all of my objects. I get performance issues for a few seconds due to probably loading the tree and stack-overflow error at x>~700 objects(mono-develop points to the insert function of the BarnesHutTree class when overflow happens). I think the tree's recursion might never end so its quickly using up all memory.

And Yes..Real Time

Im going to try out the link TiG, Thanks!

Spoiler: Code

Code (CSharp):

1. using UnityEngine;
2. using System.Collections;
3. using BHTree;
4.  
5. namespace BHTree
6. {
7.     public class BarnesHutTree
8.     {
9.    
10.         public double Theta = 0.1;
11.    
12.         public Body body;
13.         public Quad quad;
14.         public BarnesHutTree NW;
15.         public BarnesHutTree NE;
16.         public BarnesHutTree SW;
17.         public BarnesHutTree SE;
18.    
19.         public BarnesHutTree(Quad quad)
20.         {
21.             this.quad = quad;
22.             this.body = null;
23.             this.NW = null;
24.             this.NE = null;
25.             this.SW = null;
26.             this.SE = null;
27.         }
28.    
29.         public void insert(Body b)
30.         {
31.             if(this.body == null)
32.             {
33.                 body = b;
34.                 return;
35.             }
36.        
37.             if(! isExternal())
38.             {
39.                 body = body.plus(b);
40.            
41.                 putBody(b);
42.             }
43.             else
44.             {
45.                 NW = new BarnesHutTree(quad.NW());
46.                 NE = new BarnesHutTree(quad.NE());
47.                 SE = new BarnesHutTree(quad.SE());
48.                 SW = new BarnesHutTree(quad.SW());
49.            
50.                 // recursively insert both this body and Body b into the appropriate quadrant
51.                 putBody(this.body);
52.                 putBody(b);
53.                 // update the center-of-mass and total mass
54.                 body = body.plus(b);
55.             }
56.         }
57.    
58.         public void putBody(Body b)
59.         {
60.             if (b.inCheck(quad.NW()))
61.                 NW.insert(b);
62.             else if (b.inCheck(quad.NE()))
63.                 NE.insert(b);
64.             else if (b.inCheck(quad.SE()))
65.                 SE.insert(b);
66.             else if (b.inCheck(quad.SW()))
67.                 SW.insert(b);
68.         }
69.    
70.    
71.         public bool isExternal() {
72.             // a node is external iff all four children are null
73.             return (NW == null && NE == null && SW == null && SE == null);
74.         }
75.  
76.         public void updateForce(Body b)
77.         {
78.             if (body == null || b.Equals(body))
79.                 return;
80.             // if the current node is external, update net force acting on b
81.             if (isExternal())
82.                 b.addForce(body);
83.             // for internal nodes
84.             else {
85.                 // width of region represented by internal node
86.                 double s = quad.getLength();
87.                 // distance between Body b and this node's center-of-mass
88.                 double d = body.distanceTo(b);
89.                 // compare ratio (s / d) to threshold value Theta
90.                 if ((s / d) < Theta)
91.                     b.addForce(body); // b is far away
92.                 // recurse on each of current node's children
93.                 else {
94.                     NW.updateForce(b);
95.                     NE.updateForce(b);
96.                     SW.updateForce(b);
97.                     SE.updateForce(b);
98.                 }
99.             }
100.         }
101.     }
102.  
103.     public class Quad
104.     {
105.    
106.         public double xmid;
107.         public double ymid;
108.         public double length;
109.    
110.         public Quad(double xmid, double ymid, double length)
111.         {
112.             this.xmid = xmid;
113.             this.ymid = ymid;
114.             this.length = length;
115.         }
116.    
117.         public double getLength()
118.         {
119.             return length;
120.         }
121.    
122.         public bool contains(double x, double y)
123.         {
124.             double halfLen = this.length / 2.0;
125.             return (x <= this.xmid + halfLen &&
126.                     x >= this.xmid - halfLen &&
127.                     y <= this.ymid + halfLen &&
128.                     y >= this.ymid - halfLen);
129.         }
130.    
131.         public Quad NW()
132.         {
133.             double x = this.xmid - this.length / 4.0;
134.             double y = this.ymid + this.length / 4.0;
135.             double len = this.length / 2.0;
136.             Quad NW = new Quad(x, y, len);
137.             return NW;
138.         }
139.    
140.         public Quad NE()
141.         {
142.             double x = this.xmid + this.length / 4.0;
143.             double y = this.ymid + this.length / 4.0;
144.             double len = this.length / 2.0;
145.             Quad NE = new Quad(x, y, len);
146.             return NE;
147.         }
148.    
149.         public Quad SW()
150.         {
151.             double x = this.xmid - this.length / 4.0;
152.             double y = this.ymid - this.length / 4.0;
153.             double len = this.length / 2.0;
154.             Quad SW = new Quad(x, y, len);
155.             return SW;
156.         }
157.    
158.         public Quad SE()
159.         {
160.             double x = this.xmid + this.length / 4.0;
161.             double y = this.ymid - this.length / 4.0;
162.             double len = this.length / 2.0;
163.             Quad SE = new Quad(x, y, len);
164.             return SE;
165.         }
166.     }
167.  
168.     public class Body
169.     {
170.    
171.         public double G;
172.    
173.         public double rx, ry; // position
174.         public double fx, fy; // force
175.         public double mass; // mass
176.    
177.         public Body(double rx, double ry, double mass, double G)
178.         {
179.             this.rx = rx;
180.             this.ry = ry;
181.             this.mass = mass;
182.             this.G = G;
183.         }
184.    
185.         public void update(double rx, double ry)
186.         {
187.             this.rx = rx;
188.             this.ry = ry;
189.         }
190.    
191.         public double distanceTo(Body b)
192.         {
193.             double dx = rx - b.rx;
194.             double dy = ry - b.ry;
195.             return Mathf.Sqrt((float)(dx*dx + dy*dy));
196.         }
197.    
198.         public void resetForce()
199.         {
200.             fx = 0.0;
201.             fy = 0.0;
202.         }
203.    
204.         public void addForce(Body b)
205.         {
206.             if(distanceTo(b) > 0.75f)
207.             {
208.                 Body a = this;
209.                 double dx = b.rx - a.rx;
210.                 double dy = b.ry - a.ry;
211.                 double dist = Mathf.Sqrt((float)(dx*dx + dy*dy));
212.                 double F = (G * a.mass * b.mass) / (dist*dist);
213.                 a.fx += F * dx / dist;
214.                 a.fy += F * dy / dist;
215.             }
216.         }
217.    
218.         public bool inCheck(Quad q)
219.         {
220.             return q.contains(this.rx, this.ry);
221.         }
222.    
223.         public Body plus(Body b)
224.         {
225.             Body a = this;
226.             double m = a.mass + b.mass;
227.             double x = (a.rx * a.mass + b.rx * b.mass) / m;
228.             double y = (a.ry * a.mass + b.ry * b.mass) / m;
229.             return new Body(x, y, m, a.G);
230.         }
231.     }
232. }

Spoiler: More Code

Code (CSharp):

1. using UnityEngine;
2. using System.Collections;
3. using System.Collections.Generic;
4. using BHTree;
5.  
6. public class NBody : MonoBehaviour {
7.  
8.     public GameObject planetPrefab01;
9.  
10.     public double dt = 0.1;
11.     public double radius = 1000;
12.     public int N = 5;
13.     public double G = 1;
14.  
15.     private Body[] nBodies;
16.     private GameObject[] pBodies;
17.  
18.     void Start()
19.     {
20.         Initialize();
21.     }
22.  
23.     void FixedUpdate()
24.     {
25.         Quad quad = new Quad(0, 0, radius * 2);
26.         BarnesHutTree tree = new BarnesHutTree(quad);
27.         pBodies = GameObject.FindGameObjectsWithTag("Celestial Object");
28.  
29.         for(int i = 0; i < N; i++)
30.         {
31.             if(nBodies[i].inCheck(quad))
32.             {
33.                 tree.insert(nBodies[i]);
34.             }
35.         }
36.  
37.         for(int i = 0; i < N; i++)
38.         {
39.             nBodies[i].resetForce();
40.             tree.updateForce(nBodies[i]);
41.             nBodies[i].update(pBodies[i].transform.position.x, pBodies[i].transform.position.y);
42.         }
43.  
44.         for(int i = 0; i < N; i++)
45.         {
46.             Vector2 force = new Vector2((float)nBodies[i].fx, (float)nBodies[i].fy);
47.             pBodies[i].rigidbody2D.AddForce(force);
48.         }
49.     }
50.  
51.     void Initialize()
52.     {
53.         List<int> places = new List<int>();
54.         bool check = true;
55.         for(int i = 0; i < N; i++)
56.         {
57.             check = true;
58.             while(check)
59.             {
60.                 int x = Random.Range(0,(int)radius);
61.                 int y = Random.Range(0,(int)radius);
62.                 if(!places.Contains(x + y))
63.                 {
64.                     places.Add(x + y);
65.                     Instantiate(planetPrefab01, new Vector3(x, y) , Quaternion.identity);
66.                     check = false;
67.                 }
68.             }
69.         }
70.  
71.         pBodies = GameObject.FindGameObjectsWithTag("Celestial Object");
72.         nBodies = new Body[N];
73.  
74.         for(int i = 0; i < N; i++)
75.         {
76.             double px = pBodies[i].transform.position.x;
77.             double py = pBodies[i].transform.position.y;
78.             double mass = pBodies[i].rigidbody2D.mass;
79.             nBodies[i] = new Body(px, py, mass, G);
80.         }
81.     }
82. }
83.