kernel3d.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 "kernel3d.hpp"
00020 
00021 //double Kernel3d::_mindif = 1e-8;
00022 
00023 // ---------------------------------------------------------------------- 
00024 int Kernel3d::srcDOF() const
00025 {
00026   int dof = 0;
00027   switch(_kernelType) {
00028          //laplace kernels
00029   case KNL_LAP_S_U: dof = 1; break;
00030   case KNL_LAP_D_U: dof = 1; break;
00031   case KNL_LAP_I  : dof = 1; break;
00032          //stokes kernels
00033   case KNL_STK_F_U: dof = 4; break;
00034   case KNL_STK_S_U: dof = 3; break;
00035   case KNL_STK_S_P: dof = 3; break;
00036   case KNL_STK_D_U: dof = 3; break;
00037   case KNL_STK_D_P: dof = 3; break;
00038   case KNL_STK_R_U: dof = 3; break;
00039   case KNL_STK_R_P: dof = 3; break;
00040   case KNL_STK_I  : dof = 3; break;
00041   case KNL_STK_E  : dof = 3; break;
00042          //navier kernels:       //case KNL_NAV_F_U: dof = 3; break; //used for fmm
00043   case KNL_NAV_S_U: dof = 3; break;
00044   case KNL_NAV_D_U: dof = 3; break;
00045   case KNL_NAV_R_U: dof = 3; break;
00046   case KNL_NAV_I  : dof = 3; break;
00047   case KNL_NAV_E  : dof = 3; break;
00048          //others
00049   case KNL_SQRTLAP: dof = 1; break;
00050   case KNL_EXP:     dof = 1; break;
00051          //error
00052   case KNL_ERR:     dof = 0; break;
00053   }
00054   return dof;
00055 }
00056 
00057 // ---------------------------------------------------------------------- 
00058 int Kernel3d::trgDOF() const
00059 {
00060   int dof = 0;
00061   switch(_kernelType) {
00062          //laplace kernels
00063   case KNL_LAP_S_U: dof = 1; break;
00064   case KNL_LAP_D_U: dof = 1; break;
00065   case KNL_LAP_I  : dof = 1; break;
00066          //stokes kernels
00067   case KNL_STK_F_U: dof = 3; break;
00068   case KNL_STK_S_U: dof = 3; break;
00069   case KNL_STK_S_P: dof = 1; break;
00070   case KNL_STK_D_U: dof = 3; break;
00071   case KNL_STK_D_P: dof = 1; break;
00072   case KNL_STK_R_U: dof = 3; break;
00073   case KNL_STK_R_P: dof = 1; break;
00074   case KNL_STK_I  : dof = 3; break;
00075   case KNL_STK_E  : dof = 3; break;
00076          //navier kernels:       //  case KNL_NAV_F_U: dof = 3; break; //used for fmm
00077   case KNL_NAV_S_U: dof = 3; break;
00078   case KNL_NAV_D_U: dof = 3; break;
00079   case KNL_NAV_R_U: dof = 3; break;
00080   case KNL_NAV_I  : dof = 3; break;
00081   case KNL_NAV_E  : dof = 3; break;
00082          //others
00083   case KNL_SQRTLAP: dof = 1; break;
00084   case KNL_EXP    : dof = 1; break;
00085          //error
00086   case KNL_ERR:     dof = 0; break;
00087   }
00088   return dof;
00089 }
00090 
00091 // ---------------------------------------------------------------------- 
00092 bool Kernel3d::homogeneous() const
00093 {
00094   bool ret = false;
00095   switch(_kernelType) {
00096          //laplace kernels
00097   case KNL_LAP_S_U: ret = true; break;
00098          //stokes kernels
00099   case KNL_STK_F_U: ret = true; break;
00100          //navier kernels
00101   case KNL_NAV_S_U: ret = true; break;
00102          //others:
00103   case KNL_SQRTLAP: ret = true; break;
00104   case KNL_EXP    : ret = false; break;
00105   default: assert(0);
00106   }
00107   return ret;
00108 }
00109 
00110 // ---------------------------------------------------------------------- 
00111 void Kernel3d::homogeneousDeg(vector<double>& degreeVec) const
00112 {
00113   switch(_kernelType) {
00114          //laplace
00115   case KNL_LAP_S_U: degreeVec.resize(1); degreeVec[0]=1; break;
00116          //stokes
00117   case KNL_STK_F_U: degreeVec.resize(4); degreeVec[0]=1; degreeVec[1]=1; degreeVec[2]=1; degreeVec[3]=2; break;
00118          //navier
00119   case KNL_NAV_S_U: degreeVec.resize(3); degreeVec[0]=1; degreeVec[1]=1; degreeVec[2]=1; break;
00120          //others
00121   case KNL_SQRTLAP: degreeVec.resize(1); degreeVec[0]=0.5; break;
00122   case KNL_EXP    : degreeVec.resize(0); break;
00123   default: assert(0);
00124   }
00125   return;
00126 }
00127 
00128 // kernel function should be split up
00129 // ---------------------------------------------------------------------- 
00130 int Kernel3d::kernel(const DblNumMat& srcPos, const DblNumMat& srcNor, const DblNumMat& trgPos, DblNumMat& inter)
00131 {
00132   int srcDOF = this->srcDOF();
00133   int trgDOF = this->trgDOF();
00134   iA(srcPos.m()==dim() && srcNor.m()==dim() && trgPos.m()==dim());
00135   iA(srcPos.n()==srcNor.n());
00136   iA(srcPos.n()*srcDOF == inter.n());
00137   iA(trgPos.n()*trgDOF == inter.m());
00138   if(       _kernelType==KNL_LAP_S_U) {
00139          //----------------------------------------------------------------------------------
00140          //---------------------------------
00142          double OOFP = 1.0/(4.0*M_PI);
00143          for(int i=0; i<trgPos.n(); i++) {
00144                 for(int j=0; j<srcPos.n(); j++) {
00145                   double x = trgPos(0,i) - srcPos(0,j);
00146                   double y = trgPos(1,i) - srcPos(1,j);
00147                   double z = trgPos(2,i) - srcPos(2,j);
00148                   double r2 = x*x + y*y + z*z;
00149                   double r = sqrt(r2);
00150                   if(r<_mindif) {
00151                          inter(i,j) = 0;
00152                   } else {
00153                          inter(i,j) = OOFP / r;
00154                   }
00155                 }
00156          }
00157   } else if(_kernelType==KNL_LAP_D_U) {
00158          //---------------------------------
00160          double OOFP = 1.0/(4.0*M_PI);
00161          for(int i=0; i<trgPos.n(); i++)
00162                 for(int j=0; j<srcPos.n(); j++) {
00163                   double x = trgPos(0,i) - srcPos(0,j);
00164                   double y = trgPos(1,i) - srcPos(1,j);
00165                   double z = trgPos(2,i) - srcPos(2,j);
00166                   double r2 = x*x + y*y + z*z;
00167                   double r = sqrt(r2);
00168                   if(r<_mindif) {
00169                          for(int t=0;t<trgDOF;t++)
00170                                 for(int s=0;s<srcDOF;s++)
00171                                   { inter(i*trgDOF+t, j*srcDOF+s) = 0.0; }
00172                   } else {
00173                          double nx = srcNor(0,j);
00174                          double ny = srcNor(1,j);
00175                          double nz = srcNor(2,j);
00176                          double rn = x*nx + y*ny + z*nz;
00177                          double r3 = r2*r;
00178                          inter(i,j) = - OOFP / r3 * rn;
00179                   }
00180                 }
00181   } else if(_kernelType==KNL_LAP_I) {
00182          //---------------------------------
00183          for(int i=0; i<trgPos.n(); i++)
00184                 for(int j=0; j<srcPos.n(); j++)
00185                   inter(i,j) = 1;
00186   } else if(_kernelType==KNL_STK_F_U) {
00187          //----------------------------------------------------------------------------------
00188          //---------------------------------
00189          iA(_coefs.size()>=1);
00190          double mu = _coefs[0];
00191          double OOFP = 1.0/(4.0*M_PI);
00192          double OOEP = 1.0/(8.0*M_PI);
00193          double oomu = 1.0/mu;
00194          for(int i=0; i<trgPos.n(); i++)
00195                 for(int j=0; j<srcPos.n(); j++) {
00196                   double x = trgPos(0,i) - srcPos(0,j);
00197                   double y = trgPos(1,i) - srcPos(1,j);
00198                   double z = trgPos(2,i) - srcPos(2,j);
00199                   double r2 = x*x + y*y + z*z;
00200                   double r = sqrt(r2);
00201                   if(r<_mindif) {
00202                          for(int t=0;t<trgDOF;t++)
00203                                 for(int s=0;s<srcDOF;s++)
00204                                   { inter(i*trgDOF+t, j*srcDOF+s) = 0.0; }
00205                   } else {
00206                          double r3 = r2*r;
00207                          double G = oomu * OOEP / r;
00208                          double H = oomu * OOEP / r3;
00209                          double A = OOFP / r3;
00210                          int is = i*3;                   int js = j*4;
00211                          inter(is,   js)   = G + H*x*x; inter(is,   js+1) =     H*x*y; inter(is,   js+2) =     H*x*z; inter(is,   js+3) = A*x;
00212                          inter(is+1, js)   =     H*y*x; inter(is+1, js+1) = G + H*y*y; inter(is+1, js+2) =     H*y*z; inter(is+1, js+3) = A*y;
00213                          inter(is+2, js)   =     H*z*x; inter(is+2, js+1) =     H*z*y; inter(is+2, js+2) = G + H*z*z; inter(is+2, js+3) = A*z;
00214                   }
00215                 }
00216   } else if(_kernelType==KNL_STK_S_U) {
00217          //---------------------------------
00218          iA(_coefs.size()>=1);
00219          double mu = _coefs[0];
00221          double OOEP = 1.0/(8.0*M_PI);
00222          double oomu = 1.0/mu;
00223          for(int i=0; i<trgPos.n(); i++)
00224                 for(int j=0; j<srcPos.n(); j++) {
00225                   double x = trgPos(0,i) - srcPos(0,j);
00226                   double y = trgPos(1,i) - srcPos(1,j);
00227                   double z = trgPos(2,i) - srcPos(2,j);
00228                   double r2 = x*x + y*y + z*z;
00229                   double r = sqrt(r2);
00230                   if(r<_mindif) {
00231                          for(int t=0;t<trgDOF;t++)
00232                                 for(int s=0;s<srcDOF;s++)
00233                                   { inter(i*trgDOF+t, j*srcDOF+s) = 0.0; }
00234                   } else {
00235                          double r3 = r2*r;
00236                          double G = OOEP / r;
00237                          double H = OOEP / r3;
00238                          inter(i*3,   j*3)   = oomu*(G + H*x*x); inter(i*3,   j*3+1) = oomu*(    H*x*y); inter(i*3,   j*3+2) = oomu*(    H*x*z);
00239                          inter(i*3+1, j*3)   = oomu*(    H*y*x); inter(i*3+1, j*3+1) = oomu*(G + H*y*y); inter(i*3+1, j*3+2) = oomu*(    H*y*z);
00240                          inter(i*3+2, j*3)   = oomu*(    H*z*x); inter(i*3+2, j*3+1) = oomu*(    H*z*y); inter(i*3+2, j*3+2) = oomu*(G + H*z*z);
00241                   }
00242                 }
00243   } else if(_kernelType==KNL_STK_S_P) {
00244          //---------------------------------
00245          iA(_coefs.size()>=1);
00246          double mu = _coefs[0];
00248          double OOFP = 1.0/(4.0*M_PI);
00249          for(int i=0; i<trgPos.n(); i++)
00250                 for(int j=0; j<srcPos.n(); j++) {
00251                   double x = trgPos(0,i) - srcPos(0,j);
00252                   double y = trgPos(1,i) - srcPos(1,j);
00253                   double z = trgPos(2,i) - srcPos(2,j);
00254                   double r2 = x*x + y*y + z*z;
00255                   double r = sqrt(r2);
00256                   if(r<_mindif) {
00257                          for(int t=0;t<trgDOF;t++)
00258                                 for(int s=0;s<srcDOF;s++)
00259                                   { inter(i*trgDOF+t, j*srcDOF+s) = 0.0; }
00260                   } else {
00261                          double r3 = r2*r;
00262                          inter(i  ,j*3  ) = OOFP*x/r3;                   inter(i  ,j*3+1) = OOFP*y/r3;                   inter(i  ,j*3+2) = OOFP*z/r3;
00263                   }
00264                 }
00265   } else if(_kernelType==KNL_STK_D_U) {
00266          //---------------------------------
00267          iA(_coefs.size()>=1);
00268          double mu = _coefs[0];
00270          double SOEP = 6.0/(8.0*M_PI);
00271          for(int i=0; i<trgPos.n(); i++)
00272                 for(int j=0; j<srcPos.n(); j++) {
00273                   double x = trgPos(0,i) - srcPos(0,j);
00274                   double y = trgPos(1,i) - srcPos(1,j);
00275                   double z = trgPos(2,i) - srcPos(2,j);
00276                   double r2 = x*x + y*y + z*z;
00277                   double r = sqrt(r2);
00278                   if(r<_mindif) {
00279                          for(int t=0;t<trgDOF;t++)
00280                                 for(int s=0;s<srcDOF;s++)
00281                                   { inter(i*trgDOF+t, j*srcDOF+s) = 0.0; }
00282                   } else {
00283                          double nx = srcNor(0,j);
00284                          double ny = srcNor(1,j);
00285                          double nz = srcNor(2,j);
00286                          double rn = x*nx + y*ny + z*nz;
00287                          double r5 = r2*r2*r;
00288                          double C = - SOEP / r5;
00289                          inter(i*3,   j*3)   = C*rn*x*x;                                  inter(i*3,   j*3+1) = C*rn*x*y;                                 inter(i*3,   j*3+2) = C*rn*x*z;
00290                          inter(i*3+1, j*3)   = C*rn*y*x;                                  inter(i*3+1, j*3+1) = C*rn*y*y;                                 inter(i*3+1, j*3+2) = C*rn*y*z;
00291                          inter(i*3+2, j*3)   = C*rn*z*x;                                  inter(i*3+2, j*3+1) = C*rn*z*y;                                 inter(i*3+2, j*3+2) = C*rn*z*z;
00292                   }
00293                 }
00294   } else if(_kernelType==KNL_STK_D_P) {
00295          //---------------------------------
00296          iA(_coefs.size()>=1);
00297          double mu = _coefs[0];
00299          double OOTP = 1.0/(2.0*M_PI);
00300          double coef = mu*OOTP;
00301          for(int i=0; i<trgPos.n(); i++)
00302                 for(int j=0; j<srcPos.n(); j++) {
00303                   double x = trgPos(0,i) - srcPos(0,j);
00304                   double y = trgPos(1,i) - srcPos(1,j);
00305                   double z = trgPos(2,i) - srcPos(2,j);
00306                   double r2 = x*x + y*y + z*z;
00307                   double r = sqrt(r2);
00308                   if(r<_mindif) {
00309                          for(int t=0;t<trgDOF;t++)
00310                                 for(int s=0;s<srcDOF;s++)
00311                                   { inter(i*trgDOF+t, j*srcDOF+s) = 0.0; }
00312                   } else {
00313                          double nx = srcNor(0,j);
00314                          double ny = srcNor(1,j);
00315                          double nz = srcNor(2,j);
00316                          double rn = x*nx + y*ny + z*nz;
00317                          double r3 = r2*r;
00318                          double r5 = r3*r2;
00319                          int is = i;                     int js = j*3;
00320                          inter(is  ,js  ) = coef*(nx/r3 - 3*rn*x/r5);
00321                          inter(is  ,js+1) = coef*(ny/r3 - 3*rn*y/r5);
00322                          inter(is  ,js+2) = coef*(nz/r3 - 3*rn*z/r5);
00323                   }
00324                 }
00325   } else if(_kernelType==KNL_STK_R_U) {
00326          //---------------------------------
00327          iA(_coefs.size()>=1);
00328          double mu = _coefs[0];
00329          for(int i=0; i<trgPos.n(); i++)
00330                 for(int j=0; j<srcPos.n(); j++) {
00331                   double x = trgPos(0,i) - srcPos(0,j);
00332                   double y = trgPos(1,i) - srcPos(1,j);
00333                   double z = trgPos(2,i) - srcPos(2,j);
00334                   double r2 = x*x + y*y + z*z;
00335                   double r = sqrt(r2);
00336                   if(r<_mindif) {
00337                          for(int t=0;t<trgDOF;t++)
00338                                 for(int s=0;s<srcDOF;s++)
00339                                   { inter(i*trgDOF+t, j*srcDOF+s) = 0.0; }
00340                   } else {
00341                          double r3 = r2*r;
00342                          double coef = 1.0/(8.0*M_PI*mu)/r3;
00343                          inter(i*3,   j*3)   = coef*0;                   inter(i*3,   j*3+1) = coef*z;                   inter(i*3,   j*3+2) = coef*(-y);
00344                          inter(i*3+1, j*3)   = coef*(-z);                inter(i*3+1, j*3+1) = coef*0;                   inter(i*3+1, j*3+2) = coef*x;
00345                          inter(i*3+2, j*3)   = coef*y;                   inter(i*3+2, j*3+1) = coef*(-x);                inter(i*3+2, j*3+2) = coef*0;
00346                   }
00347                 }
00348   } else if(_kernelType==KNL_STK_R_P) {
00349          //---------------------------------
00350          iA(_coefs.size()>=1);
00351          double mu = _coefs[0];
00352          for(int i=0; i<trgPos.n(); i++)
00353                 for(int j=0; j<srcPos.n(); j++) {
00354                   inter(i,j*3  ) = 0;
00355                   inter(i,j*3+1) = 0;
00356                   inter(i,j*3+2) = 0;
00357                 }
00358   } else if(_kernelType==KNL_STK_I) {
00359          //---------------------------------
00360          iA(_coefs.size()>=1);
00361          double mu = _coefs[0];
00362          for(int i=0; i<trgPos.n(); i++)
00363                 for(int j=0; j<srcPos.n(); j++) {
00364                   inter(i*3,   j*3  ) = 1;               inter(i*3,   j*3+1) = 0;                inter(i*3,   j*3+2) = 0;
00365                   inter(i*3+1, j*3  ) = 0;               inter(i*3+1, j*3+1) = 1;                inter(i*3+1, j*3+2) = 0;
00366                   inter(i*3+2, j*3  ) = 0;               inter(i*3+2, j*3+1) = 0;                inter(i*3+2, j*3+2) = 1;
00367                 }
00368   } else if(_kernelType==KNL_STK_E) { //levi-civita tensor
00369          //---------------------------------
00370          iA(_coefs.size()>=1);
00371          double mu = _coefs[0];
00372          for(int i=0; i<trgPos.n(); i++)
00373                 for(int j=0; j<srcPos.n(); j++) {
00374                   double x = trgPos(0,i) - srcPos(0,j);
00375                   double y = trgPos(1,i) - srcPos(1,j);
00376                   double z = trgPos(2,i) - srcPos(2,j);
00377                   inter(i*3,   j*3  ) = 0;               inter(i*3,   j*3+1) = z;                inter(i*3,   j*3+2) = -y;
00378                   inter(i*3+1, j*3  ) = -z;      inter(i*3+1, j*3+1) = 0;                inter(i*3+1, j*3+2) = x;
00379                   inter(i*3+2, j*3  ) = y;               inter(i*3+2, j*3+1) = -x;               inter(i*3+2, j*3+2) = 0;
00380                 }
00381   } else if(_kernelType==KNL_NAV_S_U) {
00382          //----------------------------------------------------------------------------------
00383          //---------------------------------
00384          iA(_coefs.size()>=2);
00385          double mu = _coefs[0];
00386          double ve = _coefs[1];
00387          double sc1 = (3.0-4.0*ve)/(16.0*M_PI*(1.0-ve));
00388          double sc2 = 1.0/(16.0*M_PI*(1.0-ve));
00389          double oomu = 1.0/mu;
00390          for(int i=0; i<trgPos.n(); i++)
00391                 for(int j=0; j<srcPos.n(); j++) {
00392                   double x = trgPos(0,i) - srcPos(0,j);
00393                   double y = trgPos(1,i) - srcPos(1,j);
00394                   double z = trgPos(2,i) - srcPos(2,j);
00395                   double r2 = x*x + y*y + z*z;
00396                   double r = sqrt(r2);
00397                   if(r<_mindif) {
00398                          for(int t=0;t<trgDOF;t++)
00399                                 for(int s=0;s<srcDOF;s++)
00400                                   { inter(i*trgDOF+t, j*srcDOF+s) = 0.0; }
00401                   } else {
00402                          double r3 = r2*r;
00403                          double G = sc1 / r;
00404                          double H = sc2 / r3;
00405                          inter(i*3,   j*3)   = oomu*(G + H*x*x);  inter(i*3,   j*3+1) = oomu*(    H*x*y);  inter(i*3,   j*3+2) = oomu*(    H*x*z);
00406                          inter(i*3+1, j*3)   = oomu*(    H*y*x);  inter(i*3+1, j*3+1) = oomu*(G + H*y*y);  inter(i*3+1, j*3+2) = oomu*(    H*y*z);
00407                          inter(i*3+2, j*3)   = oomu*(    H*z*x);  inter(i*3+2, j*3+1) = oomu*(    H*z*y);  inter(i*3+2, j*3+2) = oomu*(G + H*z*z);
00408                   }
00409                 }
00410   } else if(_kernelType==KNL_NAV_D_U) {
00411          //---------------------------------
00412          iA(_coefs.size()>=2);
00413          double mu = _coefs[0];
00414          double ve = _coefs[1];
00415          double dc1 = -(1-2.0*ve)/(8.0*M_PI*(1.0-ve));
00416          double dc2 =  (1-2.0*ve)/(8.0*M_PI*(1.0-ve));
00417          double dc3 = -3.0/(8.0*M_PI*(1.0-ve));
00418          for(int i=0; i<trgPos.n(); i++)
00419                 for(int j=0; j<srcPos.n(); j++) {
00420                   double x = trgPos(0,i) - srcPos(0,j);
00421                   double y = trgPos(1,i) - srcPos(1,j);
00422                   double z = trgPos(2,i) - srcPos(2,j);
00423                   double r2 = x*x + y*y + z*z;
00424                   double r = sqrt(r2);
00425                   if(r<_mindif) {
00426                          for(int t=0;t<trgDOF;t++)
00427                                 for(int s=0;s<srcDOF;s++)
00428                                   { inter(i*trgDOF+t, j*srcDOF+s) = 0.0; }
00429                   } else {
00430                          double nx = srcNor(0,j);
00431                          double ny = srcNor(1,j);
00432                          double nz = srcNor(2,j);
00433                          double rn = x*nx + y*ny + z*nz;
00434                          double r3 = r2*r;
00435                          double r5 = r3*r2;
00436                          double A = dc1 / r3;                    double B = dc2 / r3;                    double C = dc3 / r5;
00437                          double rx = x;                  double ry = y;                  double rz = z;
00438                          //&&&&&&&
00439                          inter(i*3,   j*3)   = A*(rn+nx*rx) + B*(rx*nx) + C*rn*rx*rx;
00440                          inter(i*3,   j*3+1) = A*(   nx*ry) + B*(rx*ny) + C*rn*rx*ry;
00441                          inter(i*3,   j*3+2) = A*(   nx*rz) + B*(rx*nz) + C*rn*rx*rz;
00442                          inter(i*3+1, j*3)   = A*(   ny*rx) + B*(ry*nx) + C*rn*ry*rx;
00443                          inter(i*3+1, j*3+1) = A*(rn+ny*ry) + B*(ry*ny) + C*rn*ry*ry;
00444                          inter(i*3+1, j*3+2) = A*(   ny*rz) + B*(ry*nz) + C*rn*ry*rz;
00445                          inter(i*3+2, j*3)   = A*(   nz*rx) + B*(rz*nx) + C*rn*rz*rx;
00446                          inter(i*3+2, j*3+1) = A*(   nz*ry) + B*(rz*ny) + C*rn*rz*ry;
00447                          inter(i*3+2, j*3+2) = A*(rn+nz*rz) + B*(rz*nz) + C*rn*rz*rz;
00448                   }
00449                 }
00450   } else if(_kernelType==KNL_NAV_R_U) {
00451          //---------------------------------
00452          iA(_coefs.size()>=2);
00453          double mu = _coefs[0];
00454          double ve = _coefs[1];
00455          for(int i=0; i<trgPos.n(); i++)
00456                 for(int j=0; j<srcPos.n(); j++) {
00457                   double x = trgPos(0,i) - srcPos(0,j);
00458                   double y = trgPos(1,i) - srcPos(1,j);
00459                   double z = trgPos(2,i) - srcPos(2,j);
00460                   double r2 = x*x + y*y + z*z;
00461                   double r = sqrt(r2);
00462                   if(r<_mindif) {
00463                          for(int t=0;t<trgDOF;t++)
00464                                 for(int s=0;s<srcDOF;s++)
00465                                   { inter(i*trgDOF+t, j*srcDOF+s) = 0.0; }
00466                   } else {
00467                          double r3 = r2*r;
00468                          double coef = 1.0/(8.0*M_PI*mu)/r3;
00469                          inter(i*3,   j*3)   = coef*0;                   inter(i*3,   j*3+1) = coef*z;                   inter(i*3,   j*3+2) = coef*(-y);
00470                          inter(i*3+1, j*3)   = coef*(-z);                inter(i*3+1, j*3+1) = coef*0;                   inter(i*3+1, j*3+2) = coef*x;
00471                          inter(i*3+2, j*3)   = coef*y;                   inter(i*3+2, j*3+1) = coef*(-x);                inter(i*3+2, j*3+2) = coef*0;
00472                   }
00473                 }
00474   } else if(_kernelType==KNL_NAV_I) {
00475          //---------------------------------
00476          iA(_coefs.size()>=2);
00477          double mu = _coefs[0];
00478          double ve = _coefs[1];
00479          for(int i=0; i<trgPos.n(); i++)
00480                 for(int j=0; j<srcPos.n(); j++) {
00481                   inter(i*3,   j*3  ) = 1;               inter(i*3,   j*3+1) = 0;                inter(i*3,   j*3+2) = 0;
00482                   inter(i*3+1, j*3  ) = 0;               inter(i*3+1, j*3+1) = 1;                inter(i*3+1, j*3+2) = 0;
00483                   inter(i*3+2, j*3  ) = 0;               inter(i*3+2, j*3+1) = 0;                inter(i*3+2, j*3+2) = 1;
00484                 }
00485   } else if(_kernelType==KNL_NAV_E) {
00486          //---------------------------------
00487          iA(_coefs.size()>=2);
00488          double mu = _coefs[0];
00489          double ve = _coefs[1];
00490          for(int i=0; i<trgPos.n(); i++)
00491                 for(int j=0; j<srcPos.n(); j++) {
00492                   double x = trgPos(0,i) - srcPos(0,j);
00493                   double y = trgPos(1,i) - srcPos(1,j);
00494                   double z = trgPos(2,i) - srcPos(2,j);
00495                   inter(i*3,   j*3  ) = 0;               inter(i*3,   j*3+1) = z;                inter(i*3,   j*3+2) = -y;
00496                   inter(i*3+1, j*3  ) = -z;      inter(i*3+1, j*3+1) = 0;                inter(i*3+1, j*3+2) = x;
00497                   inter(i*3+2, j*3  ) = y;               inter(i*3+2, j*3+1) = -x;               inter(i*3+2, j*3+2) = 0;
00498                 }
00499   } else if(_kernelType==KNL_SQRTLAP) {
00500          //----------------------------------------------------------------------------------
00501          //---------------------------------
00502          for(int i=0; i<trgPos.n(); i++) {
00503                 for(int j=0; j<srcPos.n(); j++) {
00504                   double x = trgPos(0,i) - srcPos(0,j);
00505                   double y = trgPos(1,i) - srcPos(1,j);
00506                   double z = trgPos(2,i) - srcPos(2,j);
00507                   double r2 = x*x + y*y + z*z;
00508                   double r = sqrt(r2);
00509                   if(r<_mindif) {
00510                          inter(i,j) = 0;
00511                   } else {
00512                          inter(i,j) = 1 / sqrt(r);
00513                   }
00514                 }
00515          }
00516   } else if(_kernelType==KNL_EXP) {
00517          //---------------------------------
00518          for(int i=0; i<trgPos.n(); i++) {
00519                 for(int j=0; j<srcPos.n(); j++) {
00520                   double x = trgPos(0,i) - srcPos(0,j);
00521                   double y = trgPos(1,i) - srcPos(1,j);
00522                   double z = trgPos(2,i) - srcPos(2,j);
00523                   double r2 = x*x + y*y + z*z;
00524                   double r = sqrt(r2);
00525                   if(r<_mindif) {
00526                          inter(i,j) = 0;
00527                   } else {
00528                          inter(i,j) = exp(-r2);
00529                   }
00530                 }
00531          }
00532   } else if(_kernelType==KNL_ERR) {
00533          //----------------------------------------------------------------------------------
00534          //---------------------------------
00535          iA(0);
00536   }
00537   return 0;
00538 }

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