orsa_interaction_tree.cc

Go to the documentation of this file.
00001 /* 
00002    ORSA - Orbit Reconstruction, Simulation and Analysis
00003    Copyright (C) 2002-2004 Pasquale Tricarico
00004    
00005    This program is free software; you can redistribute it and/or
00006    modify it under the terms of the GNU General Public License
00007    as published by the Free Software Foundation; either version 2
00008    of the License, or (at your option) any later version.
00009    
00010    As a special exception, Pasquale Tricarico gives permission to
00011    link this program with Qt commercial edition, and distribute the
00012    resulting executable, without including the source code for the Qt
00013    commercial edition in the source distribution.
00014    
00015    This program is distributed in the hope that it will be useful,
00016    but WITHOUT ANY WARRANTY; without even the implied warranty of
00017    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00018    GNU General Public License for more details.
00019    
00020    You should have received a copy of the GNU General Public License
00021    along with this program; if not, write to the Free Software
00022    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
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       // parent = 0;
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     // TreeNode *parent;
00075     // cubic domains... *** o is included, o+l is included ***
00076     Vector o; // cube origin (minimum x,y,z)
00077     double l; // cube side
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     // monopole
00105     
00106     const double l2 = d.LengthSquared();
00107     
00108     if (d.IsZero()) {
00109       cerr << "*** Warning: two objects in the same position! (" << l2 << ")" << endl;
00110       // continue;
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     // quadrupole
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)); // works?
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     // printf("inside: (%g,%g,%g)      origin: (%g,%g,%g) l: %g\n",p.x,p.y,p.z,o.x,o.y,o.z,l);
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)); // works?
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     // zero bodies
00263     if (b.begin() == b.end()) return;
00264   
00265     // one body
00266     if (++b.begin() == b.end()) return;
00267     
00268     // compute domain for root
00269     if (root) { 
00270       depth = 0;
00271       
00272       Vector p; // (maximum x,y,z)
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       // printf("total bodies: %i\n",total_bodies);
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     // a slightly bigger root domain (NEEDED!)...
00301     //  if (root) { // only for root?
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     // check
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     // clear child list
00326     child.clear();
00327     
00328     // build eight nodes
00329     {
00330       TreeNode n;
00331       n.l = l/2.0;
00332       n.depth = depth+1;
00333       // n.parent = this;
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     // fill new nodes with bodies...
00346     {
00347       list<TreeNode>::iterator c_it;
00348       list<Body>::iterator b_it = b.begin();
00349       while (b_it != b.end()) {
00350         /* 
00351            printf("body position: (%g,%g,%g)\n",
00352            b_it->position().x,
00353            b_it->position().y,
00354            b_it->position().z);
00355         */
00356         
00357         c_it = child.begin();
00358         while (c_it != child.end()) {
00359           if (c_it->inside_domain(b_it->position())) {
00360             // printf("--- inside domain!!! ---\n");
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     // remove empty childs
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     // refine mesh
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(); // Note: this *may* be O(N)
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     // root_node.print();
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               // cerr << "num_domain: " << num_domain << "  num_direct: " << num_direct << "  ratio: " << (1.0*num_domain)/(1.0*num_direct) << endl;
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                 // cerr << "num_domain: " << num_domain << "  num_direct: " << num_direct << "  ratio: " << (1.0*num_domain)/(1.0*num_direct) << endl;
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 { // not leaf
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     // to be done...
00562     return 0.0;
00563   }
00564   
00565 } // namespace orsa

Generated on Fri Nov 3 20:37:42 2006 for liborsa by  doxygen 1.4.7