fmm3d.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 "common/vecmatop.hpp"
00019 #include "fmm3d.hpp"
00020 
00021 using std::istringstream;
00022 
00023 FMM3d::FMM3d(const string& p):
00024   KnlMat3d(p), _center(0,0,0), _rootLevel(0),
00025   _np(6), _let(NULL), _matmgnt(NULL)
00026 {
00027 }
00028 
00029 FMM3d::~FMM3d()
00030 {
00031   if(_let!=NULL)         delete _let;
00032 }
00033 
00034 // ---------------------------------------------------------------------- 
00035 int FMM3d::SrcEqu2TrgChk_dgemv(const DblNumMat& srcPos, const DblNumMat& srcNor, const DblNumMat& trgPos, const DblNumVec& srcDen, DblNumVec& trgVal)
00036 {
00037   int TMAX = 1024;
00038   if(trgPos.n()<=TMAX) {
00039          int M = trgPos.n() * _knl.trgDOF();
00040          int N = srcPos.n() * _knl.srcDOF();
00041          DblNumMat tmp(M,N);
00042          iC( _knl.kernel(srcPos, srcNor, trgPos, tmp) );
00043          iC( dgemv(1.0, tmp, srcDen, 1.0, trgVal) );
00044   } else {
00045          int RUNS = (trgPos.n()-1) / TMAX + 1;
00046          for(int r=0; r<RUNS; r++) {
00047                 int stt = r*TMAX;
00048                 int end = min((r+1)*TMAX, trgPos.n());
00049                 int num = end-stt;
00050                 int M = num * _knl.trgDOF();
00051                 int N = srcPos.n() * _knl.srcDOF();
00052                 DblNumMat tps(dim(), num, false, trgPos.data() + stt*dim() );
00053                 DblNumVec tvl(num*_knl.trgDOF(), false, trgVal.data() + stt*_knl.trgDOF());
00054                 DblNumMat tmp(M,N);
00055                 iC( _knl.kernel(srcPos, srcNor, tps, tmp) );
00056                 iC( dgemv(1.0, tmp, srcDen, 1.0, tvl) );
00057          }
00058   }
00059   return (0);
00060 }
00061 // ---------------------------------------------------------------------- 
00062 int FMM3d::SrcEqu2UpwChk_dgemv(const DblNumMat& srcPos, const DblNumMat& srcNor, Point3 trgCtr, double trgRad, const DblNumVec& srcDen, DblNumVec& trgVal)
00063 {
00064   DblNumMat trgPos; iC( _matmgnt->localPos(UC, trgCtr, trgRad, trgPos) );
00065   int M = trgPos.n() * _knl.trgDOF();
00066   int N = srcPos.n() * _knl.srcDOF();
00067   DblNumMat tmp(M,N);
00068   iC( _knl.kernel(srcPos, srcNor, trgPos, tmp) );
00069   iC( dgemv(1.0, tmp, srcDen, 1.0, trgVal) );
00070   return (0);
00071 }
00072 // ---------------------------------------------------------------------- 
00073 int FMM3d::SrcEqu2DwnChk_dgemv(const DblNumMat& srcPos, const DblNumMat& srcNor, Point3 trgCtr, double trgRad, const DblNumVec& srcDen, DblNumVec& trgVal)
00074 {
00075   DblNumMat trgPos; iC( _matmgnt->localPos(DC, trgCtr, trgRad, trgPos) );
00076   int M = trgPos.n() * _knl.trgDOF();
00077   int N = srcPos.n() * _knl.srcDOF();
00078   DblNumMat tmp(M,N);
00079   iC( _knl.kernel(srcPos, srcNor, trgPos, tmp) );
00080   iC( dgemv(1.0, tmp, srcDen, 1.0, trgVal) );
00081   return 0;
00082 }
00083 // ---------------------------------------------------------------------- 
00084 int FMM3d::DwnEqu2TrgChk_dgemv(Point3 srcCtr, double srcRad, const DblNumMat& trgPos, const DblNumVec& srcDen, DblNumVec& trgVal)
00085 {
00086   int TMAX = 1024;
00087   if(trgPos.n()<=TMAX) {
00088          DblNumMat srcPos; iC( _matmgnt->localPos(DE, srcCtr, srcRad, srcPos) );
00089          int M = trgPos.n() * _knl_mm.trgDOF();
00090          int N = srcPos.n() * _knl_mm.srcDOF();
00091          DblNumMat tmp(M,N);
00092          iC( _knl_mm.kernel(srcPos, srcPos, trgPos, tmp) );
00093          iC( dgemv(1.0, tmp, srcDen, 1.0, trgVal) );
00094   } else {
00095          DblNumMat srcPos; iC( _matmgnt->localPos(DE, srcCtr, srcRad, srcPos) );
00096          int RUNS = (trgPos.n()-1) / TMAX + 1;
00097          for(int r=0; r<RUNS; r++) {
00098                 int stt = r*TMAX;
00099                 int end = min((r+1)*TMAX, trgPos.n());
00100                 int num = end-stt;
00101                 int M = num * _knl_mm.trgDOF();
00102                 int N = srcPos.n() * _knl_mm.srcDOF();
00103                 DblNumMat tps(dim(), num, false, trgPos.data() + stt*dim());
00104                 DblNumVec tvl(num*_knl_mm.trgDOF(), false, trgVal.data() + stt*_knl_mm.trgDOF());
00105                 DblNumMat tmp(M, N);
00106                 iC( _knl_mm.kernel(srcPos, srcPos, tps, tmp) );
00107                 iC( dgemv(1.0, tmp, srcDen, 1.0, tvl) );
00108          }
00109   }
00110   return 0;
00111 }
00112 
00113 // ---------------------------------------------------------------------- 
00114 int FMM3d::UpwEqu2TrgChk_dgemv(Point3 srcCtr, double srcRad, const DblNumMat& trgPos, const DblNumVec& srcDen, DblNumVec& trgVal)
00115 {
00116   int TMAX = 1024;
00117   if(trgPos.n()<=TMAX) {
00118          DblNumMat srcPos; iC( _matmgnt->localPos(UE, srcCtr, srcRad, srcPos) );
00119          int M = trgPos.n() * _knl_mm.trgDOF();
00120          int N = srcPos.n() * _knl_mm.srcDOF();
00121          DblNumMat tmp(M,N);
00122          iC( _knl_mm.kernel(srcPos, srcPos, trgPos, tmp) );
00123          iC( dgemv(1.0, tmp, srcDen, 1.0, trgVal) );
00124   } else {
00125          DblNumMat srcPos; iC( _matmgnt->localPos(UE, srcCtr, srcRad, srcPos) );
00126          int RUNS = (trgPos.n()-1) / TMAX + 1;
00127          for(int r=0; r<RUNS; r++) {
00128                 int stt = r*TMAX;
00129                 int end = min((r+1)*TMAX, trgPos.n());
00130                 int num = end-stt;
00131                 int M = num * _knl_mm.trgDOF();
00132                 int N = srcPos.n() * _knl_mm.srcDOF();
00133                 DblNumMat tps(dim(), num, false, trgPos.data() + stt*dim());
00134                 DblNumVec tvl(num*_knl_mm.trgDOF(), false, trgVal.data() + stt*_knl_mm.trgDOF());
00135                 DblNumMat tmp(M,N);
00136                 iC( _knl_mm.kernel(srcPos, srcPos, tps, tmp) );
00137                 iC( dgemv(1.0, tmp, srcDen, 1.0, tvl) );
00138          }
00139   }
00140   return 0;
00141 }
00142 
00143 // ---------------------------------------------------------------------- 
00144 DblNumMat FMM3d::srcExaPos(int gNodeIdx) 
00145 {
00146   Let3d::Node& node=_let->node(gNodeIdx);
00147   int beg = node.srcExaBeg();
00148   int num = node.srcExaNum();
00149   return DblNumMat(dim(), num, false, _srcExaPos.data()+beg*dim());
00150 }
00151 DblNumMat FMM3d::srcExaNor(int gNodeIdx)
00152 {
00153   Let3d::Node& node=_let->node(gNodeIdx);
00154   int beg = node.srcExaBeg();
00155   int num = node.srcExaNum();
00156   return DblNumMat(dim(), num, false, _srcExaNor.data()+beg*dim());
00157 }
00158 DblNumVec FMM3d::srcExaDen(int gNodeIdx)
00159 {
00160   Let3d::Node& node=_let->node(gNodeIdx);
00161   int beg = node.srcExaBeg();
00162   int num = node.srcExaNum();
00163   return DblNumVec(srcDOF()*num, false, _srcExaDen.data()+beg*srcDOF());
00164 }
00165 DblNumVec FMM3d::srcUpwEquDen(int gNodeIdx)
00166 {
00167   Let3d::Node& node=_let->node(gNodeIdx);
00168   int idx = node.srcNodeIdx();
00169   return DblNumVec(datSze(UE), false, _srcUpwEquDen.data()+idx*datSze(UE));
00170 }
00171 DblNumVec FMM3d::srcUpwChkVal(int gNodeIdx)
00172 {
00173   Let3d::Node& node=_let->node(gNodeIdx);
00174   int idx = node.srcNodeIdx();
00175   return DblNumVec(datSze(UC), false, _srcUpwChkVal.data()+idx*datSze(UC));
00176 }
00177 // ---------------------------------------------------------------------- 
00178 DblNumMat FMM3d::trgExaPos(int gNodeIdx)
00179 {
00180   Let3d::Node& node=_let->node(gNodeIdx);
00181   int beg = node.trgExaBeg();
00182   int num = node.trgExaNum();
00183   return DblNumMat(dim(), num, false, _trgExaPos.data()+beg*dim());
00184 }
00185 DblNumVec FMM3d::trgExaVal(int gNodeIdx)
00186 {
00187   Let3d::Node& node=_let->node(gNodeIdx);
00188   int beg = node.trgExaBeg();
00189   int num = node.trgExaNum();
00190   return DblNumVec(trgDOF()*num, false, _trgExaVal.data()+beg*trgDOF());
00191 }
00192 DblNumVec FMM3d::trgDwnEquDen(int gNodeIdx)
00193 {
00194   Let3d::Node& node=_let->node(gNodeIdx);
00195   int idx = node.trgNodeIdx();
00196   return DblNumVec(datSze(DE), false, _trgDwnEquDen.data()+idx*datSze(DE));
00197 }
00198 DblNumVec FMM3d::trgDwnChkVal(int gNodeIdx)
00199 {
00200   Let3d::Node& node=_let->node(gNodeIdx);
00201   int idx = node.trgNodeIdx();
00202   return DblNumVec(datSze(DC), false, _trgDwnChkVal.data()+idx*datSze(DC));
00203 }
00204 
00205 

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