00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #include "orsa_interaction.h"
00026 #include "orsa_secure_math.h"
00027 #include "orsa_universe.h"
00028
00029 #include <iostream>
00030 #include <list>
00031 #include <stack>
00032 #include <map>
00033
00034 using namespace std;
00035
00036 namespace orsa {
00037
00038 class TreeNode {
00039 public:
00040 TreeNode() {
00041 reset();
00042 }
00043
00044 void reset() {
00045 bool_node_mass_computed = false;
00046 bool_node_quadrupole_computed = false;
00047 bool_node_center_of_mass_computed = false;
00048
00049 depth = 0;
00050 }
00051
00052 bool inside_domain(const Vector p) const;
00053
00054 bool is_leaf() const {
00055 return (child.empty() && (!b.empty()));
00056 };
00057
00058 double node_mass() const;
00059
00060 double * node_quadrupole() const;
00061
00062 Vector node_center_of_mass() const;
00063
00064 void BuildMesh(const bool root=false);
00065
00066 public:
00067 void print() const;
00068
00069 public:
00070 list<Body> b;
00071
00072 public:
00073 list<TreeNode> child;
00074
00075
00076 Vector o;
00077 double l;
00078 unsigned int depth;
00079
00080 private:
00081 mutable double _node_mass;
00082 mutable bool bool_node_mass_computed;
00083
00084 mutable double _node_quadrupole[3][3];
00085 mutable bool bool_node_quadrupole_computed;
00086
00087 mutable Vector _node_center_of_mass;
00088 mutable bool bool_node_center_of_mass_computed;
00089 };
00090
00091 double delta_function(const unsigned int i, const unsigned int j) {
00092 if (i==j) return 1.0;
00093 return 0.0;
00094 }
00095
00096 Vector ComputeAcceleration(const list<Body>::const_iterator body_it, const list<TreeNode>::const_iterator node_domain_it, const bool compute_quadrupole=true) {
00097
00098 Vector a;
00099
00100 if (node_domain_it->node_mass()==0.0) return a;
00101
00102 Vector d = node_domain_it->node_center_of_mass() - body_it->position();
00103
00104
00105
00106 const double l2 = d.LengthSquared();
00107
00108 if (d.IsZero()) {
00109 cerr << "*** Warning: two objects in the same position! (" << l2 << ")" << endl;
00110
00111 return a;
00112 }
00113
00114 a += d * secure_pow(l2,-1.5) * node_domain_it->node_mass();
00115
00116 if (!compute_quadrupole) {
00117 return a;
00118 }
00119
00120
00121
00122 double x[3];
00123
00124 x[0] = d.x;
00125 x[1] = d.y;
00126 x[2] = d.z;
00127
00128 double coefficient = 0.0;
00129 unsigned int i,j;
00130 double c_node_quadrupole[3][3];
00131 memcpy((void*)c_node_quadrupole, (const void*)node_domain_it->node_quadrupole(), 3*3*sizeof(double));
00132 for (i=0;i<3;++i) {
00133 for (j=0;j<3;++j) {
00134 coefficient += c_node_quadrupole[i][j]*x[i]*x[j];
00135 }
00136 }
00137
00138 a += d * secure_pow(l2,-3.0) * coefficient;
00139
00140 return a;
00141 }
00142
00143 bool TreeNode::inside_domain(const Vector p) const {
00144
00145
00146
00147 if (p.x < o.x) return false;
00148 if (p.y < o.y) return false;
00149 if (p.z < o.z) return false;
00150
00151 if (p.x > o.x+l) return false;
00152 if (p.y > o.y+l) return false;
00153 if (p.z > o.z+l) return false;
00154
00155 return true;
00156 }
00157
00158 double TreeNode::node_mass() const {
00159 if (bool_node_mass_computed) return _node_mass;
00160 _node_mass = 0.0;
00161 {
00162 list<TreeNode>::const_iterator c_it=child.begin();
00163 while (c_it!=child.end()) {
00164 _node_mass += c_it->node_mass();
00165 ++c_it;
00166 }
00167 }
00168 {
00169 list<Body>::const_iterator b_it=b.begin();
00170 while (b_it!=b.end()) {
00171 _node_mass += b_it->mass();
00172 ++b_it;
00173 }
00174 }
00175 bool_node_mass_computed = true;
00176 return _node_mass;
00177 }
00178
00179 double * TreeNode::node_quadrupole() const {
00180 if (bool_node_quadrupole_computed) return (&_node_quadrupole[0][0]);
00181 unsigned int i,j;
00182 for (i=0;i<3;++i) {
00183 for (j=0;j<3;++j) {
00184 _node_quadrupole[i][j] = 0.0;
00185 }
00186 }
00187 double x[3];
00188 double l_sq;
00189 double c_node_quadrupole[3][3];
00190 Vector vec;
00191 {
00192 list<TreeNode>::const_iterator c_it=child.begin();
00193 while (c_it!=child.end()) {
00194 vec = c_it->node_center_of_mass() - node_center_of_mass();
00195
00196 x[0] = vec.x;
00197 x[1] = vec.y;
00198 x[2] = vec.z;
00199
00200 l_sq = vec.LengthSquared();
00201
00202 memcpy((void*)c_node_quadrupole, (const void*)c_it->node_quadrupole(), 3*3*sizeof(double));
00203
00204 for (i=0;i<3;++i) {
00205 for (j=0;j<3;++j) {
00206 _node_quadrupole[i][j] += c_it->node_mass()*(3.0*x[i]*x[j]-l_sq*delta_function(i,j)) + c_node_quadrupole[i][j];
00207 }
00208 }
00209 ++c_it;
00210 }
00211 }
00212 {
00213 list<Body>::const_iterator b_it=b.begin();
00214 while (b_it!=b.end()) {
00215 vec = b_it->position() - node_center_of_mass();
00216
00217 x[0] = vec.x;
00218 x[1] = vec.y;
00219 x[2] = vec.z;
00220
00221 l_sq = vec.LengthSquared();
00222
00223 for (i=0;i<3;++i) {
00224 for (j=0;j<3;++j) {
00225 _node_quadrupole[i][j] += b_it->mass()*(3.0*x[i]*x[j]-l_sq*delta_function(i,j));
00226 }
00227 }
00228 ++b_it;
00229 }
00230 }
00231 bool_node_mass_computed = true;
00232 return (&_node_quadrupole[0][0]);
00233 }
00234
00235 Vector TreeNode::node_center_of_mass() const {
00236 if (bool_node_center_of_mass_computed) return _node_center_of_mass;
00237 Vector vec_sum;
00238 double mass_sum=0;
00239 {
00240 list<TreeNode>::const_iterator c_it=child.begin();
00241 while (c_it!=child.end()) {
00242 vec_sum += c_it->node_mass()*c_it->node_center_of_mass();
00243 mass_sum += c_it->node_mass();
00244 ++c_it;
00245 }
00246 }
00247 {
00248 list<Body>::const_iterator b_it=b.begin();
00249 while (b_it!=b.end()) {
00250 vec_sum += b_it->mass()*b_it->position();
00251 mass_sum += b_it->mass();
00252 ++b_it;
00253 }
00254 }
00255 _node_center_of_mass = vec_sum/mass_sum;
00256 bool_node_center_of_mass_computed = true;
00257 return _node_center_of_mass;
00258 }
00259
00260 void TreeNode::BuildMesh(const bool root) {
00261
00262
00263 if (b.begin() == b.end()) return;
00264
00265
00266 if (++b.begin() == b.end()) return;
00267
00268
00269 if (root) {
00270 depth = 0;
00271
00272 Vector p;
00273 o = p = b.begin()->position();
00274 Vector r;
00275 list<Body>::iterator b_it = b.begin();
00276 ++b_it;
00277 unsigned int total_bodies=1;
00278 while (b_it != b.end()) {
00279 r = b_it->position();
00280
00281 if (r.x<o.x) o.x = r.x;
00282 if (r.y<o.y) o.y = r.y;
00283 if (r.z<o.z) o.z = r.z;
00284
00285 if (r.x>p.x) p.x = r.x;
00286 if (r.y>p.y) p.y = r.y;
00287 if (r.z>p.z) p.z = r.z;
00288
00289 ++total_bodies;
00290 ++b_it;
00291 }
00292
00293
00294
00295 l = (p.x-o.x);
00296 if ((p.y-o.y)>l) l = (p.y-o.y);
00297 if ((p.z-o.z)>l) l = (p.z-o.z);
00298 }
00299
00300
00301
00302 {
00303 const double d=0.01*l;
00304
00305 o.x -= d;
00306 o.y -= d;
00307 o.z -= d;
00308
00309 l += 2.0*d;
00310 }
00311
00312
00313 {
00314 list<Body>::iterator b_it = b.begin();
00315 while (b_it != b.end()) {
00316 if (inside_domain(b_it->position())) {
00317
00318 } else {
00319 printf("WARNING! One body outside domain...\n");
00320 }
00321 ++b_it;
00322 }
00323 }
00324
00325
00326 child.clear();
00327
00328
00329 {
00330 TreeNode n;
00331 n.l = l/2.0;
00332 n.depth = depth+1;
00333
00334
00335 n.o.Set(o.x ,o.y ,o.z ); child.push_back(n);
00336 n.o.Set(o.x ,o.y ,o.z+n.l); child.push_back(n);
00337 n.o.Set(o.x ,o.y+n.l,o.z ); child.push_back(n);
00338 n.o.Set(o.x ,o.y+n.l,o.z+n.l); child.push_back(n);
00339 n.o.Set(o.x+n.l,o.y ,o.z ); child.push_back(n);
00340 n.o.Set(o.x+n.l,o.y ,o.z+n.l); child.push_back(n);
00341 n.o.Set(o.x+n.l,o.y+n.l,o.z ); child.push_back(n);
00342 n.o.Set(o.x+n.l,o.y+n.l,o.z+n.l); child.push_back(n);
00343 }
00344
00345
00346 {
00347 list<TreeNode>::iterator c_it;
00348 list<Body>::iterator b_it = b.begin();
00349 while (b_it != b.end()) {
00350
00351
00352
00353
00354
00355
00356
00357 c_it = child.begin();
00358 while (c_it != child.end()) {
00359 if (c_it->inside_domain(b_it->position())) {
00360
00361 c_it->b.push_back(*b_it);
00362
00363 b.erase(b_it);
00364 b_it = b.begin();
00365 if (b_it == b.end()) break;
00366
00367 c_it = child.begin();
00368 --c_it;
00369 }
00370 ++c_it;
00371 }
00372 ++b_it;
00373 }
00374 }
00375
00376
00377 {
00378 list<TreeNode>::iterator c_it = child.begin();
00379 while (c_it != child.end()) {
00380 if (c_it->b.empty()) {
00381 c_it = child.erase(c_it);
00382 } else {
00383 ++c_it;
00384 }
00385 }
00386 }
00387
00388
00389 {
00390 list<TreeNode>::iterator c_it = child.begin();
00391 while (c_it != child.end()) {
00392 c_it->BuildMesh();
00393 ++c_it;
00394 }
00395 }
00396 }
00397
00398 void TreeNode::print() const {
00399 unsigned int bodies=0;
00400 list<Body>::const_iterator b_it = b.begin();
00401 while (b_it != b.end()) {
00402 ++bodies;
00403 ++b_it;
00404 }
00405
00406 unsigned int childs=child.size();
00407
00408 printf("node --- depth: %i childs: %i mass: %g cube side: %g origin: (%g,%g,%g) bodies: %i\n",depth,childs,node_mass(),l,o.x,o.y,o.z,bodies);
00409
00410 list<TreeNode>::const_iterator it=child.begin();
00411 while (it!=child.end()) {
00412 it->print();
00413 ++it;
00414 }
00415 }
00416
00417 GravitationalTree::GravitationalTree() : Interaction() {
00418 g = GetG();
00419 theta = 0.7;
00420 }
00421
00422 GravitationalTree::GravitationalTree(const GravitationalTree &) : Interaction() {
00423 g = GetG();
00424 theta = 0.7;
00425 }
00426
00427 Interaction * GravitationalTree::clone() const {
00428 return new GravitationalTree(*this);
00429 }
00430
00431 void GravitationalTree::Acceleration(const Frame &f, vector<Vector> &a) {
00432
00433 if (f.size() < 2) return;
00434
00435 a.resize(f.size());
00436
00437 {
00438 unsigned int k=0;
00439 while (k<a.size()) {
00440 a[k].Set(0.0,0.0,0.0);
00441 ++k;
00442 }
00443 }
00444
00445 map <unsigned int, unsigned int> frame_map;
00446 {
00447 unsigned int k=0;
00448 while (k<f.size()) {
00449 frame_map[f[k].BodyId()] = k;
00450 ++k;
00451 }
00452 }
00453
00454 TreeNode root_node;
00455
00456 unsigned int k=0;
00457 while (k<f.size()) {
00458 root_node.b.push_back(f[k]);
00459 ++k;
00460 }
00461
00462 root_node.BuildMesh(true);
00463
00464
00465
00466 list<TreeNode>::const_iterator node_body_it, node_domain_it;
00467 list<Body>::const_iterator body_it;
00468
00469 stack<list<TreeNode>::const_iterator> stk_body, stk_domain;
00470
00471 double angle;
00472
00473 unsigned int num_direct=0, num_domain=0;
00474
00475 node_body_it = root_node.child.begin();
00476 while (node_body_it != root_node.child.end()) {
00477
00478 if (node_body_it->is_leaf()) {
00479
00480 body_it = node_body_it->b.begin();
00481 while (body_it != node_body_it->b.end()) {
00482
00483 node_domain_it = root_node.child.begin();
00484 while (node_domain_it != root_node.child.end()) {
00485
00486 angle = (node_domain_it->l)/(node_domain_it->node_center_of_mass()-body_it->position()).Length();
00487
00488 if (angle < theta) {
00489
00490 ++num_domain;
00491
00492
00493 a[frame_map[body_it->BodyId()]] += ComputeAcceleration(body_it,node_domain_it);
00494
00495 ++node_domain_it;
00496
00497 } else if (node_domain_it->is_leaf()) {
00498
00499 if (body_it->BodyId() != node_domain_it->b.begin()->BodyId()) {
00500 a[frame_map[body_it->BodyId()]] += ComputeAcceleration(body_it,node_domain_it);
00501 ++num_direct;
00502
00503 }
00504
00505 ++node_domain_it;
00506
00507 } else {
00508
00509 stk_domain.push(node_domain_it);
00510 node_domain_it = node_domain_it->child.begin();
00511
00512 }
00513
00514 while (stk_domain.size()) {
00515 if (node_domain_it == stk_domain.top()->child.end()) {
00516 node_domain_it = stk_domain.top();
00517 ++node_domain_it;
00518 stk_domain.pop();
00519 } else {
00520 break;
00521 }
00522 }
00523
00524 }
00525
00526 ++body_it;
00527 }
00528
00529 ++node_body_it;
00530
00531 } else {
00532
00533 stk_body.push(node_body_it);
00534 node_body_it = node_body_it->child.begin();
00535
00536 }
00537
00538 while (stk_body.size()) {
00539 if (node_body_it == stk_body.top()->child.end()) {
00540 node_body_it = stk_body.top();
00541 ++node_body_it;
00542 stk_body.pop();
00543 } else {
00544 break;
00545 }
00546 }
00547
00548 }
00549
00550 {
00551 unsigned int k=0;
00552 while (k<a.size()) {
00553 a[k] *= g;
00554 ++k;
00555 }
00556 }
00557
00558 }
00559
00560 double GravitationalTree::PotentialEnergy(const Frame&) {
00561
00562 return 0.0;
00563 }
00564
00565 }