let3d.cpp

Go to the documentation of this file.
00001 /* Kernel Independent Fast Multipole Method
00002    Copyright (C) 2004 Lexing Ying, New York University
00003 
00004 This program is free software; you can redistribute it and/or modify
00005 it under the terms of the GNU General Public License as published by
00006 the Free Software Foundation; either version 2, or (at your option)
00007 any later version.
00008 
00009 This program is distributed in the hope that it will be useful, but WITHOUT
00010 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00011 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00012 for more details.
00013 
00014 You should have received a copy of the GNU General Public License
00015 along with this program; see the file COPYING.  If not, write to the Free
00016 Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
00017 02111-1307, USA.  */
00018 #include "let3d.hpp"
00019 
00020 using std::min;
00021 using std::max;
00022 using std::set;
00023 using std::queue;
00024 using std::ofstream;
00025 using std::cerr;
00026 using std::endl;
00027 using std::istringstream;
00028 
00029 //-----------------------------------------------------------------------
00030 Let3d::Let3d(const string& p):
00031   ComObject(p), _srcPos(NULL), _trgPos(NULL), _center(0.0), _rootLevel(0),
00032   _ptsMax(150), _maxLevel(10)
00033 {
00034 }
00035 
00036 Let3d::~Let3d()
00037 {
00038 }
00039 
00040 // ---------------------------------------------------------------------- 
00041 int Let3d::setup(map<string,string>& optionsMap)
00042 {
00043   //------------
00044   map<string,string>::iterator mapindex;
00045   mapindex = optionsMap.find("-" + prefix() + "ptsMax"); iA(mapindex!=optionsMap.end());
00046   { istringstream ss((*mapindex).second); ss>>_ptsMax; }
00047   mapindex = optionsMap.find("-" + prefix() + "maxLevel"); iA(mapindex!=optionsMap.end());
00048   { istringstream ss((*mapindex).second); ss>>_maxLevel; }
00049   //------------
00050   iC( srcData() );
00051   iC( trgData() );
00052   return 0;
00053 }
00054 
00055 // ---------------------------------------------------------------------- 
00056 int Let3d::srcData()
00057 {
00058   //-----------------------------------------
00060   DblNumMat& pos = *(_srcPos);  iA( pos.m()==dim() );
00061   
00062   vector<Node>& nodeVec = this->nodeVec(); nodeVec.clear();
00063   vector< vector<int> > vecIdxs;  
00064   //local src number, the number of src point in each box
00065   vector<int> lclSrcNumVec; //glb 
00066   
00067   nodeVec.push_back( Node(-1,-1, Index3(0,0,0), 0) );
00068   vecIdxs.push_back( vector<int>() );
00069   vector<int>& curVecIdxs = vecIdxs[0];
00070   Point3 bbmin(center()-Point3(radius()));
00071   Point3 bbmax(center()+Point3(radius()));
00072   for(int k=0; k<pos.n(); k++) {
00073          Point3 tmp(pos.clmdata(k));     
00074          iA(tmp>=bbmin && tmp<=bbmax);//LEXING: IMPORANT
00075          curVecIdxs.push_back(k);
00076   }
00077   lclSrcNumVec.push_back( curVecIdxs.size() );
00078   
00079   int level = 0;
00080   int arrBeg = 0;
00081   int arrEnd = 1;
00082   int arrCnt = 0;
00083   while(arrBeg < arrEnd) {
00084          //1.
00085          arrCnt = arrEnd;
00086          for(int k=arrBeg; k < arrEnd; k++) {
00087                 //---
00088                 if( lclSrcNumVec[k]>ptsMax() && level<maxLevel()-1 ) {
00089                   nodeVec[k].child() = arrCnt;
00090                   arrCnt = arrCnt + pow2(dim());
00091                   //children's ess                
00092                   for(int a=0; a<2; a++) {
00093                          for(int b=0; b<2; b++) {
00094                                 for(int c=0; c<2; c++) {
00095                                   nodeVec.push_back( Node(k,-1, 2*nodeVec[k].path2Node()+Index3(a,b,c), nodeVec[k].depth()+1) ); //par, chd
00096                                   vecIdxs.push_back( vector<int>() );
00097                                   lclSrcNumVec.push_back( 0 );
00098                                 }
00099                          }
00100                   }
00101                   //children's vector of indices
00102                   Point3 centerCurNode( center(k) ); //get center of current node
00103                   for(vector<int>::iterator vecIdxsIt=vecIdxs[k].begin(); vecIdxsIt!=vecIdxs[k].end(); vecIdxsIt++) {
00104                          Point3 tmp(pos.clmdata(*vecIdxsIt));
00105                          Index3 idx;
00106                          for(int j=0; j<dim(); j++){
00107                                 idx(j) = (tmp(j) >= centerCurNode(j));
00108                          }
00109                          int chdGNodeIdx = child(k, idx);
00110                          vecIdxs[chdGNodeIdx].push_back(*vecIdxsIt);
00111                   }
00112                   vecIdxs[k].clear(); //VERY IMPORTANT
00113                   //children's lsm                
00114                   for(int a=0; a<2; a++) {
00115                          for(int b=0; b<2; b++) {
00116                                 for(int c=0; c<2; c++) {
00117                                   int chdGNodeIdx = child( k, Index3(a,b,c) );
00118                                   lclSrcNumVec[chdGNodeIdx] = vecIdxs[chdGNodeIdx].size();
00119                                 }
00120                          }
00121                   }
00122                 }
00123          }
00124          level++;
00125          arrBeg = arrEnd;
00126          arrEnd = arrCnt;
00127   }
00128   _level = level; //SET LEVEL
00129 
00130   //ordering of the boxes, in top-down or bottom-up fashion
00131   vector<int> orderBoxesVec     ;  iC( dwnOrderCollect(orderBoxesVec    ) );
00132   //set other parts of essvec
00133   int cnt = 0;
00134   int sum = 0;
00135   for(int i=0; i < orderBoxesVec        .size(); i++) {
00136          int gNodeIdx = orderBoxesVec   [i];
00137          if(lclSrcNumVec[gNodeIdx]>0) {
00138                 nodeVec[gNodeIdx].tag() = nodeVec[gNodeIdx].tag() | LET_SRCNODE;
00139                 nodeVec[gNodeIdx].srcNodeIdx() = cnt;
00140                 cnt++;
00141                 if(nodeVec[gNodeIdx].child()==-1) {
00142                   nodeVec[gNodeIdx].srcExaBeg() = sum;
00143                   nodeVec[gNodeIdx].srcExaNum() = lclSrcNumVec[gNodeIdx];
00144                   sum += lclSrcNumVec[gNodeIdx];
00145                   nodeVec[gNodeIdx].srcOwnVecIdxs() = vecIdxs[gNodeIdx];
00146                 }
00147          }
00148   }
00149   _srcNodeCnt = cnt;  _srcExaCnt = sum; //SET S cnts
00150   return 0;
00151 }
00152 
00153 // ---------------------------------------------------------------------- 
00154 int Let3d::trgData()
00155 {
00156   //-----------------------------------------
00157   //edata
00158   DblNumMat& pos = *(_trgPos);  iA( pos.m()==dim() );
00159   
00160   vector<Node>& nodeVec = this->nodeVec();
00161   vector< vector<int> > vecIdxs; vecIdxs.resize(nodeVec.size() );
00162   vector<int> lclSrcNumVec;           lclSrcNumVec.resize(nodeVec.size(), 0);
00163   
00164   vector<int>& curVecIdxs = vecIdxs[0];
00165   Point3 bbmin(center()-Point3(radius()));
00166   Point3 bbmax(center()+Point3(radius()));  //cerr<<" bbx "<<bbmin<<" "<<bbmax<<endl;
00167   for(int k=0; k < pos.n(); k++) {
00168          Point3 tmp(pos.clmdata(k));
00169          iA(tmp>=bbmin && tmp<=bbmax);   //LEXING: IMPORTANT
00170          curVecIdxs.push_back(k);
00171   }
00172   lclSrcNumVec[0] = curVecIdxs.size();
00173   
00174   vector<int> orderBoxesVec     ;  iC( dwnOrderCollect(orderBoxesVec    ) );
00175   for(int i=0; i < orderBoxesVec        .size(); i++) {
00176          int gNodeIdx = orderBoxesVec[i];
00177          Node& curNode = nodeVec[gNodeIdx];
00178          vector<int>& curVecIdxs = vecIdxs[gNodeIdx];    
00179          if(curNode.child()!=-1) { //not terminal
00180                 //children's vecIdxs
00181                 Point3 curCenter( center(gNodeIdx) );
00182                 for(vector<int>::iterator curVecIdxsIt=curVecIdxs.begin();curVecIdxsIt !=curVecIdxs.end(); curVecIdxsIt++) {
00183                   Point3 tmp(pos.clmdata(*curVecIdxsIt));
00184                   Index3 idx;
00185                   for(int j=0; j<dim(); j++)
00186                          idx(j) = (tmp(j)>=curCenter(j));
00187                   int chdGNodeIdx = child(gNodeIdx, idx);
00188                   vector<int>& chdVecIdxs = vecIdxs[chdGNodeIdx];
00189                   chdVecIdxs.push_back(*curVecIdxsIt);
00190                 }
00191                 curVecIdxs.clear(); //VERY IMPORTANT
00192                 //children's lsm                //              for(int ord=0; ord<8; ord++) {
00193                 for(int a=0; a<2; a++) {
00194                   for(int b=0; b<2; b++) {
00195                          for(int c=0; c<2; c++) {
00196                                 int chdGNodeIdx = child(gNodeIdx, Index3(a,b,c));
00197                                 lclSrcNumVec[chdGNodeIdx] = vecIdxs[chdGNodeIdx].size();
00198                          }
00199                   }
00200                 }
00201          }
00202   }
00203   //set EVTR
00204   int cnt = 0;
00205   int sum = 0;
00206   for(int i=0; i<orderBoxesVec  .size(); i++) {
00207          int gNodeIdx = orderBoxesVec   [i];
00208          if(lclSrcNumVec[gNodeIdx]>0) { //evtr node
00209                 nodeVec[gNodeIdx].tag() = nodeVec[gNodeIdx].tag() | LET_TRGNODE;
00210                 nodeVec[gNodeIdx].trgNodeIdx() = cnt;
00211                 cnt ++;
00212                 if(nodeVec[gNodeIdx].child()==-1) { //terminal
00213                   nodeVec[gNodeIdx].trgExaBeg() = sum;
00214                   nodeVec[gNodeIdx].trgExaNum() = lclSrcNumVec[gNodeIdx];
00215                   sum += lclSrcNumVec[gNodeIdx];
00216                   nodeVec[gNodeIdx].trgOwnVecIdxs() = vecIdxs[gNodeIdx];
00217                 }
00218          }
00219   }
00220   _trgNodeCnt = cnt;  _trgExaCnt = sum;
00221   
00222   //set USER
00223   for(int i=0; i < orderBoxesVec        .size(); i++) {
00224          int gNodeIdx = orderBoxesVec   [i];
00225          if(nodeVec[gNodeIdx].tag() & LET_TRGNODE) { //a evtr           
00226                 iC( calgnext(gNodeIdx) );
00227          }
00228   }
00229   //cerr<<usndecnt<<" "<<usextcnt<<endl;
00230   return 0;
00231 }
00232 
00233 // ---------------------------------------------------------------------- 
00234 int Let3d::calgnext(int gNodeIdx)
00235 {
00236   vector<Node>&  nodeVec = this->nodeVec();
00237   
00238   set<int> Uset, Vset, Wset, Xset;
00239   int curGNodeIdx = gNodeIdx;
00240   if(root(curGNodeIdx)==false) {
00241          int parGNodeIdx = parent(curGNodeIdx);
00242          
00243          Index3 minIdx(0);
00244          Index3 maxIdx(pow2(depth(curGNodeIdx)));
00245 
00246          for(int i=-2; i<4; i++) {
00247                 for(int j=-2; j<4; j++) {
00248                   for(int k=-2; k<4; k++) {
00249                          Index3 tryPath( 2*path2Node(parGNodeIdx) + Index3(i,j,k) );
00250                          if(tryPath >= minIdx && tryPath <  maxIdx && tryPath != path2Node(curGNodeIdx)) {      
00251                                 int resGNodeIdx = findgnt(depth(curGNodeIdx), tryPath);
00252                                 bool adj = adjacent(resGNodeIdx, curGNodeIdx);
00253                                 if( depth(resGNodeIdx) < depth(curGNodeIdx) ) {
00254                                   if(adj){
00255                                          if(terminal(curGNodeIdx)){
00256                                                 Uset.insert(resGNodeIdx);
00257                                          }
00258                                          else { ; }
00259                                   }
00260                                   else {
00261                                          Xset.insert(resGNodeIdx);
00262                                   }
00263                                 }
00264                                 if( depth(resGNodeIdx)==depth(curGNodeIdx) ) {
00265                                   if(!adj) {
00266                                          Index3 bb(path2Node(resGNodeIdx)-path2Node(curGNodeIdx));
00267                                          assert( bb.linfty()<=3 );
00268                                          Vset.insert(resGNodeIdx);
00269                                   }
00270                                   else {
00271                                          if(terminal(curGNodeIdx)) {
00272                                                 queue<int> rest;
00273                                                 rest.push(resGNodeIdx);
00274                                                 while(rest.empty()==false) {
00275                                                   int fntGNodeIdx = rest.front(); rest.pop();                                    //int fntgNodeIdx = fntgnt.gNodeIdx();
00276                                                   if(adjacent(fntGNodeIdx, curGNodeIdx)==false) {
00277                                                          Wset.insert( fntGNodeIdx );
00278                                                   }
00279                                                   else {
00280                                                          if(terminal(fntGNodeIdx)) {
00281                                                                 Uset.insert(fntGNodeIdx);
00282                                                          }
00283                                                          else { 
00284                                                                 for(int a=0; a<2; a++) {
00285                                                                   for(int b=0; b<2; b++) {
00286                                                                          for(int c=0; c<2; c++) {
00287                                                                                 rest.push( child(fntGNodeIdx, Index3(a,b,c)) );
00288                                                                          }
00289                                                                   }
00290                                                                 }
00291                                                          }
00292                                                   }
00293                                                 }
00294                                          }
00295                                   }
00296                                 }
00297                          }
00298                   }
00299                 }
00300          }
00301   }
00302   if(terminal(curGNodeIdx))
00303          Uset.insert(curGNodeIdx);
00304   
00305   for(set<int>::iterator si=Uset.begin(); si!=Uset.end(); si++)
00306          if(nodeVec[*si].tag() & LET_SRCNODE)           nodeVec[gNodeIdx].Unodes().push_back(*si);
00307   for(set<int>::iterator si=Vset.begin(); si!=Vset.end(); si++)
00308          if(nodeVec[*si].tag() & LET_SRCNODE)           nodeVec[gNodeIdx].Vnodes().push_back(*si);
00309   for(set<int>::iterator si=Wset.begin(); si!=Wset.end(); si++)
00310          if(nodeVec[*si].tag() & LET_SRCNODE)           nodeVec[gNodeIdx].Wnodes().push_back(*si);
00311   for(set<int>::iterator si=Xset.begin(); si!=Xset.end(); si++)
00312          if(nodeVec[*si].tag() & LET_SRCNODE)           nodeVec[gNodeIdx].Xnodes().push_back(*si);
00313   
00314   return (0);
00315 }
00316 // ---------------------------------------------------------------------- 
00317 int Let3d::dwnOrderCollect(vector<int>& orderBoxesVec   )
00318 {
00319   orderBoxesVec.clear();
00320   for(int i=0; i < nodeVec().size(); i++)
00321          orderBoxesVec.push_back(i);
00322   iA(orderBoxesVec.size()==nodeVec().size());
00323   return 0;
00324 }
00325 // ---------------------------------------------------------------------- 
00326 #undef __FUNCT__
00327 #define __FUNCT__ "Let3d::upwardOrderCollect"
00328 int Let3d::upwOrderCollect(vector<int>& orderBoxesVec   )
00329 {
00330   orderBoxesVec.clear();
00331   for(int i=nodeVec().size()-1; i>=0; i--)
00332          orderBoxesVec.push_back(i);
00333   iA(orderBoxesVec.size()==nodeVec().size());
00334   return 0;
00335 }
00336 // ---------------------------------------------------------------------- 
00337 int Let3d::child(int gNodeIdx, const Index3& idx)
00338 {
00339   assert(idx>=Index3(0) && idx<Index3(2));
00340   return node(gNodeIdx).child() + (idx(0)*4+idx(1)*2+idx(2));
00341 }
00342 Point3 Let3d::center(int gNodeIdx) //center of a node
00343 {
00344   Point3 ll( center() - Point3(radius()) );
00345   int tmp = pow2(depth(gNodeIdx));
00346   Index3 pathLcl(path2Node(gNodeIdx));
00347   Point3 res;
00348   for(int d=0; d<dim(); d++) {
00349          res(d) = ll(d) + (2*radius()) * (pathLcl(d)+0.5) / double(tmp);
00350   }
00351   return res;
00352 }
00353 double Let3d::radius(int gNodeIdx) //radius of a node
00354 {
00355   return radius()/double(pow2(depth(gNodeIdx)));
00356 }
00357 // ---------------------------------------------------------------------- 
00358 int Let3d::findgnt(int wntdepth, const Index3& wntpath)
00359 {
00360   int cur = 0;  //cerr<<"GOOD "<<path(cur)<<"     ";
00361   Index3 leftpath(wntpath);
00362   while(depth(cur)<wntdepth && terminal(cur)==false) {
00363          int dif = wntdepth-depth(cur);
00364          int tmp = pow2(dif-1);
00365          Index3 choice( leftpath/tmp );
00366          leftpath -= choice*tmp;
00367          cur = child(cur, choice);       //cur = child(cur, IdxIter(0,2,true,choice) );  //cerr<<path(cur)<<"["<<choice<<" "<<tmp<<"]"<<"     ";
00368   }  //cerr<<endl;
00369   return cur;
00370 }
00371 // ---------------------------------------------------------------------- 
00372 bool Let3d::adjacent(int me, int yo)
00373 {
00374   int md = max(depth(me),depth(yo));
00375   Index3 one(1);
00376   Index3 mecenter(  (2*path2Node(me)+one) * pow2(md - depth(me))  );
00377   Index3 yocenter(  (2*path2Node(yo)+one) * pow2(md - depth(yo))  );
00378   int meradius = pow2(md - depth(me));
00379   int yoradius = pow2(md - depth(yo));
00380   Index3 dif( abs(mecenter-yocenter) );
00381   int radius  = meradius + yoradius;
00382   return
00383          ( dif <= Index3(radius) ) && //not too far
00384          ( dif.linfty() == radius ); //at least one edge touch
00385 }
00386 
00387 // ---------------------------------------------------------------------- 
00388 int Let3d::print()
00389 {
00390   cerr<<_nodeVec.size()<<" "<<_level<<endl;
00391   cerr<<_srcNodeCnt <<" "<< _srcExaCnt << endl;
00392   cerr<<_trgNodeCnt <<" "<< _trgExaCnt << endl;
00393 
00394   for(int i=0; i<_nodeVec.size(); i++) {
00395          cerr<<node(i).parent()<<" "<<node(i).child()<<endl;
00396   }
00397   
00398   for(int i=0; i<_nodeVec.size(); i++) {
00399          //cerr<<node(i).Unodes().size()<<" "<<node(i).Vnodes().size()<<" "<<node(i).Wnodes().size()<<" "<<node(i).Xnodes().size()<<" "<<endl;
00400          Node& n = node(i);
00401          cerr<<n.srcNodeIdx()<<" "<<n.srcExaBeg()<<" "<<n.srcExaNum()<<" "
00402                   <<n.trgNodeIdx()<<" "<<n.trgExaBeg()<<" "<<n.trgExaNum()<<endl;
00403   }
00404   return 0;
00405 }

Generated on Sun Dec 4 19:24:39 2005 for fmm3d by  doxygen 1.4.5