int FMM3d::evaluate const DblNumVec &  srcDen,
DblNumVec &  trgVal
[virtual]
 

Evaluate - in fmm3d_eval.cpp

Implements KnlMat3d.

Definition at line 25 of file fmm3d_eval.cpp.

References _let, _matmgnt, _srcExaDen, _srcUpwChkVal, _srcUpwEquDen, _trgDwnChkVal, _trgDwnEquDen, _trgExaVal, LET_SRCNODE, Let3d::node(), MatMgnt3d::report(), KnlMat3d::srcDOF(), srcExaDen(), Let3d::tag(), Let3d::terminal(), KnlMat3d::trgDOF(), and Let3d::upwOrderCollect().

00026 {
00027   _matmgnt->report();
00028   //-----------------------------------
00029   iA(srcDen.m()==srcDOF()*(*_srcPos).n());  iA(trgVal.m()==trgDOF()*(*_trgPos).n());
00030   
00031   //cerr<<"fmm src and trg numbers "<<pglbnum(_srcPos)<<" "<<pglbnum(_trgPos)<<endl;
00032   int srcDOF = this->srcDOF();
00033   int trgDOF = this->trgDOF();
00034   
00035   //1. zero out Vecs
00036   setvalue(trgVal, 0.0);
00037   
00038   setvalue(_srcExaDen, 0.0);
00039   setvalue(_srcUpwEquDen, 0.0);
00040   setvalue(_srcUpwChkVal, 0.0);
00041   
00042   setvalue(_trgExaVal, 0.0);
00043   setvalue(_trgDwnEquDen, 0.0);
00044   setvalue(_trgDwnChkVal, 0.0);
00045   clock_t ck0, ck1;
00046   //CLOCKING;
00047   ck0 = clock();
00048   vector<int> ordVec; iC( _let->upwOrderCollect(ordVec) ); //BOTTOM UP
00049 
00050   //2. for cbtr, load ExaDen
00051   ck0 = clock();
00052   for(int i=0; i<ordVec.size(); i++) {
00053          int gNodeIdx = ordVec[i];
00054          if(_let->tag(gNodeIdx) & LET_SRCNODE) {
00055                 if(_let->terminal(gNodeIdx)==true) {
00056                   DblNumVec srcExaDen(this->srcExaDen(gNodeIdx));
00057                   vector<int>& curVecIdxs = _let->node(gNodeIdx).srcOwnVecIdxs();
00058                   for(int k=0; k<curVecIdxs.size(); k++) {
00059                          int poff = curVecIdxs[k];
00060                          for(int d=0; d<srcDOF; d++) {
00061                                 srcExaDen(k*srcDOF+d) = srcDen(poff*srcDOF+d);
00062                          }
00063                   }
00064                 }
00065          }
00066   }
00067   //ck1 = clock();  cerr<<"load  "<<double(ck1-ck0)/CLOCKS_PER_SEC<<endl;  
00068   
00069   //3. up computation
00070   ck0 = clock();
00071   for(int i=0; i<ordVec.size(); i++) {
00072          int gNodeIdx = ordVec[i];
00073          if(_let->tag(gNodeIdx) & LET_SRCNODE) {                //GNTra gnt = _let->gNodeIdx2gnt(gNodeIdx);
00074                 //if(_let->depth(gNodeIdx)>=2) {
00075                 if(_let->depth(gNodeIdx)>=0) {
00076                   DblNumVec srcUpwChkValgNodeIdx(srcUpwChkVal(gNodeIdx));
00077                   DblNumVec srcUpwEquDengNodeIdx(srcUpwEquDen(gNodeIdx));
00078                   if(_let->terminal(gNodeIdx)==true) {
00079                          //S2M - Source -> Multipole Exapnsion
00080                          iC( SrcEqu2UpwChk_dgemv(srcExaPos(gNodeIdx), srcExaNor(gNodeIdx), _let->center(gNodeIdx), _let->radius(gNodeIdx), srcExaDen(gNodeIdx), srcUpwChkValgNodeIdx) );
00081                   } else {
00082                          //M2M - Multipole -> Multipole
00083                          for(int a=0; a<2; a++) {
00084                                 for(int b=0; b<2; b++) {
00085                                   for(int c=0; c<2; c++) {
00086                                          Index3 idx(a,b,c);
00087                                          int chi = _let->child(gNodeIdx, idx);
00088                                          if(_let->tag(chi) & LET_SRCNODE) {
00089                                                 iC( _matmgnt->UpwEqu2UpwChk_dgemv(_let->depth(chi)+_rootLevel, idx, srcUpwEquDen(chi), srcUpwChkValgNodeIdx) );
00090                                          }
00091                                   }
00092                                 }
00093                          }
00094                   }
00095                   //M2M - Multipole -> Multipole
00096                   iC( _matmgnt->UpwChk2UpwEqu_dgemv(_let->depth(gNodeIdx)+_rootLevel, srcUpwChkValgNodeIdx, srcUpwEquDengNodeIdx) );
00097                 }
00098          }
00099   }
00100   ordVec.clear();  iC( _let->dwnOrderCollect(ordVec) );
00101   //U - list contribution calculation
00102   ck0 = clock();
00103   for(int i=0; i<ordVec.size(); i++) {
00104          int gNodeIdx = ordVec[i];
00105          if(_let->tag(gNodeIdx) & LET_TRGNODE) { //evaluator
00106                 if( _let->terminal(gNodeIdx)==true ) { //terminal       
00107                   Let3d::Node& curNode = _let->node(gNodeIdx);
00108                   DblNumVec trgExaValgNodeIdx(trgExaVal(gNodeIdx));
00109                   DblNumMat trgExaPosgNodeIdx(trgExaPos(gNodeIdx));
00110                   for(vector<int>::iterator vi=curNode.Unodes().begin(); vi!=curNode.Unodes().end(); vi++) {
00111                          //S2T - source -> target
00112                          iC( SrcEqu2TrgChk_dgemv(srcExaPos(*vi), srcExaNor(*vi), trgExaPosgNodeIdx, srcExaDen(*vi), trgExaValgNodeIdx) );
00113                   }
00114                 }
00115          }
00116   }
00117   //V - list contribution calculation
00118   ck0 = clock();
00119   for(int i=0; i<ordVec.size(); i++) {
00120          int gNodeIdx = ordVec[i];
00121          if( _let->tag(gNodeIdx) & LET_TRGNODE) { //eValuator           //GNTra gnt = _let->gNodeIdx2gnt(gNodeIdx);
00122                 Point3 gNodeIdxCtr(_let->center(gNodeIdx));
00123                 double D = 2.0 * _let->radius(gNodeIdx);
00124                 DblNumVec trgDwnChkVal(this->trgDwnChkVal(gNodeIdx));
00125                 Let3d::Node& curNode = _let->node(gNodeIdx);
00126                 for(vector<int>::iterator vi=curNode.Vnodes().begin(); vi!=curNode.Vnodes().end(); vi++) {
00127                   Point3 viCtr(_let->center(*vi));
00128                   Index3 idx;
00129                   for(int d=0; d<dim(); d++){
00130                          idx(d) = int(round( (viCtr[d]-gNodeIdxCtr[d])/D ));
00131                   }
00132                   Node& srcPtr = node(*vi);
00133                   Node& trgPtr = node(gNodeIdx);
00134                   if(srcPtr.VotCnt()==0) {
00135                          srcPtr.effDen().resize( _matmgnt->effDatSze(UE) );                      setvalue(srcPtr.effDen(), 0.0);//1. resize effDen
00136                          iC( _matmgnt->plnDen2EffDen(_let->depth(gNodeIdx)+_rootLevel, srcUpwEquDen(*vi),  srcPtr.effDen()) );                   //2. transform from upeDen to effDen
00137                   }
00138                   if(trgPtr.VinCnt()==0) {
00139                          trgPtr.effVal().resize( _matmgnt->effDatSze(DC) );                      setvalue(trgPtr.effVal(), 0.0);                         //1. resize effVal
00140                   }
00141                   //M2L - multipole -> local
00142                   iC( _matmgnt->UpwEqu2DwnChk_dgemv(_let->depth(gNodeIdx)+_rootLevel, idx, srcPtr.effDen(), trgPtr.effVal()) );
00143                   
00144                   srcPtr.VotCnt()++;
00145                   trgPtr.VinCnt()++;
00146                   if(srcPtr.VotCnt()==srcPtr.VotNum()) {
00147                          srcPtr.effDen().resize(0);                      //1. resize effDen to 0
00148                          srcPtr.VotCnt()=0;
00149                   }
00150                   if(trgPtr.VinCnt()==trgPtr.VinNum()) {
00151                          iC( _matmgnt->effVal2PlnVal(_let->depth(gNodeIdx)+_rootLevel, trgPtr.effVal(), trgDwnChkVal) );                         //1. transform from effVal to dncVal
00152                          trgPtr.effVal().resize(0); //2. resize effVal to 0
00153                          trgPtr.VinCnt()=0;
00154                   }
00155                 }
00156          }
00157   }
00158   //W - list contrubtion calculation
00159   ck0 = clock();
00160   for(int i=0; i<ordVec.size(); i++) {
00161          int gNodeIdx = ordVec[i];
00162          if( _let->tag(gNodeIdx) & LET_TRGNODE ) {
00163                 if( _let->terminal(gNodeIdx)==true ) {
00164                   DblNumVec trgExaVal_gNodeIdx(this->trgExaVal(gNodeIdx));
00165                   Let3d::Node& curNode = _let->node(gNodeIdx);
00166                   for(vector<int>::iterator vi=curNode.Wnodes().begin(); vi!=curNode.Wnodes().end(); vi++) {
00167                          if(_let->terminal(*vi) && _let->node(*vi).srcExaNum()*srcDOF<_matmgnt->plnDatSze(UE)) { //use Exa instead
00168                                 //S2T - source -> target
00169                                 iC( SrcEqu2TrgChk_dgemv(srcExaPos(*vi), srcExaNor(*vi), trgExaPos(gNodeIdx), srcExaDen(*vi), trgExaVal_gNodeIdx) );
00170                          } else {
00171                                 //M2T - multipole -> target
00172                                 int vni = *vi;                          
00173                                 iC( UpwEqu2TrgChk_dgemv(_let->center(vni), _let->radius(vni), trgExaPos(gNodeIdx), srcUpwEquDen(*vi), trgExaVal_gNodeIdx) );
00174                          }
00175                   }
00176                 }
00177          }
00178   }
00179   //X - list contrubtion calculation
00180   ck0 = clock();
00181   for(int i=0; i<ordVec.size(); i++) {
00182          int gNodeIdx = ordVec[i];
00183          if( _let->tag(gNodeIdx) & LET_TRGNODE) {       
00184                 Let3d::Node& curNode = _let->node(gNodeIdx);
00185                 DblNumVec trgExaVal_gNodeIdx(trgExaVal(gNodeIdx));
00186                 DblNumVec trgDwnChkVal_gNodeIdx(trgDwnChkVal(gNodeIdx));
00187                 for(vector<int>::iterator vi=curNode.Xnodes().begin(); vi!=curNode.Xnodes().end(); vi++) {
00188                   if(_let->terminal(gNodeIdx) && _let->node(gNodeIdx).trgExaNum()*trgDOF<_matmgnt->plnDatSze(DC)) { //use Exa instead
00189                          iC( SrcEqu2TrgChk_dgemv(srcExaPos(*vi), srcExaNor(*vi), trgExaPos(gNodeIdx), srcExaDen(*vi), trgExaVal_gNodeIdx) );
00190                   } else {
00191                          //S2L - source -> local
00192                          iC( SrcEqu2DwnChk_dgemv(srcExaPos(*vi), srcExaNor(*vi), _let->center(gNodeIdx), _let->radius(gNodeIdx), srcExaDen(*vi), trgDwnChkVal_gNodeIdx) );
00193                   }
00194                 }
00195          }
00196   }
00197   //7. combine
00198   ck0 = clock();
00199   for(int i=0; i<ordVec.size(); i++) {
00200          int gNodeIdx = ordVec[i];
00201          if( _let->tag(gNodeIdx) & LET_TRGNODE ) { //eValuator          
00202                 if(_let->depth(gNodeIdx)>=3) {
00203                   int pargNodeIdx = _let->parent(gNodeIdx);     
00204                   Index3 chdIdx( _let->path2Node(gNodeIdx)-2 * _let->path2Node(pargNodeIdx) );
00205                   //L2L - local -> local
00206                   DblNumVec trgDwnChkVal_gNodeIdx(trgDwnChkVal(gNodeIdx));
00207                   iC( _matmgnt->DwnEqu2DwnChk_dgemv(_let->depth(pargNodeIdx)+_rootLevel, chdIdx, trgDwnEquDen(pargNodeIdx), trgDwnChkVal_gNodeIdx) );
00208                 }
00209                 if(_let->depth(gNodeIdx)>=2) {
00210                   //L2L - local -> local
00211                   DblNumVec trgDwnEquDen_gNodeIdx(trgDwnEquDen(gNodeIdx));
00212                   iC( _matmgnt->DwnChk2DwnEqu_dgemv(_let->depth(gNodeIdx)+_rootLevel, trgDwnChkVal(gNodeIdx), trgDwnEquDen_gNodeIdx) );
00213                 }
00214                 if(_let->terminal(gNodeIdx)) {
00215                   //L2T - local -> target
00216                   DblNumVec trgExaVal_gNodeIdx(trgExaVal(gNodeIdx));
00217                   iC( DwnEqu2TrgChk_dgemv(_let->center(gNodeIdx), _let->radius(gNodeIdx), trgExaPos(gNodeIdx), trgDwnEquDen(gNodeIdx), trgExaVal_gNodeIdx) );
00218                 }
00219          }
00220   }
00221   
00222   //8. save trgExaVal
00223   ck0 = clock();
00224   for(int i=0; i<ordVec.size(); i++) {
00225          int gNodeIdx = ordVec[i];
00226          if( _let->tag(gNodeIdx) & LET_TRGNODE ) {
00227                 if( _let->terminal(gNodeIdx)==true ) {
00228                   DblNumVec trgExaVal(this->trgExaVal(gNodeIdx));
00229                   vector<int>& curVecIdxs = _let->node(gNodeIdx).trgOwnVecIdxs();
00230                   for(int k=0; k<curVecIdxs.size(); k++) {
00231                          int poff = curVecIdxs[k];
00232                          for(int d=0; d<trgDOF; d++) {
00233                                 trgVal(poff*trgDOF+d) = trgExaVal(k*trgDOF+d);
00234                          }
00235                   }
00236                 }
00237          }
00238   }
00239   return (0);
00240 }


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