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