kernel3d_mpi.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_mpi.hpp"
00020 #include "comobject_mpi.hpp"
00021 
00030 /* Most of this code is self-explanatory.  Based on the type of kernel being used
00031  * to do the evaluations, evaluate a multiplied context and the degrees of
00032  * freedom necessary */
00033 
00034 /* Minimum difference */
00035 double Kernel3d_MPI::_mindif = 1e-8;
00036 
00037 // ----------------------------------------------------------------------
00038 /* Depending on kernel type, set the degree of freedom for the sources
00039  * See kernel3d_mpi.hpp for kernelType descriptions */
00040 int Kernel3d_MPI::srcDOF() const
00041 {
00042   int dof = 0;
00043   switch(_kernelType) {
00044          //laplace kernels
00045   case KNL_LAP_S_U: dof = 1; break;
00046   case KNL_LAP_D_U: dof = 1; break;
00047   case KNL_LAP_I  : dof = 1; break;
00048          //stokes kernels
00049   case KNL_STK_F_U: dof = 4; break;
00050   case KNL_STK_S_U: dof = 3; break;
00051   case KNL_STK_S_P: dof = 3; break;
00052   case KNL_STK_D_U: dof = 3; break;
00053   case KNL_STK_D_P: dof = 3; break;
00054   case KNL_STK_R_U: dof = 3; break;
00055   case KNL_STK_R_P: dof = 3; break;
00056   case KNL_STK_I  : dof = 3; break;
00057   case KNL_STK_E  : dof = 3; break;
00058          //navier kernels:       //case KNL_NAV_F_U: dof = 3; break; //used for fmm
00059   case KNL_NAV_S_U: dof = 3; break;
00060   case KNL_NAV_D_U: dof = 3; break;
00061   case KNL_NAV_R_U: dof = 3; break;
00062   case KNL_NAV_I  : dof = 3; break;
00063   case KNL_NAV_E  : dof = 3; break;
00064          //others
00065          //error
00066   case KNL_ERR:     dof = 0; break;
00067   }
00068   return dof;
00069 }
00070 
00071 // ----------------------------------------------------------------------
00072 /* Depending on the kernel type, set the degree of freedom for the target values
00073  * See kernel3d_mpi.hpp for kernelType descriptions */
00074 int Kernel3d_MPI::trgDOF() const
00075 {
00076   int dof = 0;
00077   switch(_kernelType) {
00078          //laplace kernels
00079   case KNL_LAP_S_U: dof = 1; break;
00080   case KNL_LAP_D_U: dof = 1; break;
00081   case KNL_LAP_I  : dof = 1; break;
00082          //stokes kernels
00083   case KNL_STK_F_U: dof = 3; break;
00084   case KNL_STK_S_U: dof = 3; break;
00085   case KNL_STK_S_P: dof = 1; break;
00086   case KNL_STK_D_U: dof = 3; break;
00087   case KNL_STK_D_P: dof = 1; break;
00088   case KNL_STK_R_U: dof = 3; break;
00089   case KNL_STK_R_P: dof = 1; break;
00090   case KNL_STK_I  : dof = 3; break;
00091   case KNL_STK_E  : dof = 3; break;
00092          //navier kernels:       //  case KNL_NAV_F_U: dof = 3; break; //used for fmm
00093   case KNL_NAV_S_U: dof = 3; break;
00094   case KNL_NAV_D_U: dof = 3; break;
00095   case KNL_NAV_R_U: dof = 3; break;
00096   case KNL_NAV_I  : dof = 3; break;
00097   case KNL_NAV_E  : dof = 3; break;
00098          //others
00099          //error
00100   case KNL_ERR:     dof = 0; break;
00101   }
00102   return dof;
00103 }
00104 
00105 // ---------------------------------------------------------------------- 
00106 bool Kernel3d_MPI::homogeneous() const
00107 {
00108   bool ret = false;
00109   switch(_kernelType) {
00110          //laplace kernels
00111   case KNL_LAP_S_U: ret = true; break;
00112          //stokes kernels
00113   case KNL_STK_F_U: ret = true; break;
00114          //stokes kernels
00115   case KNL_NAV_S_U: ret = true; break;
00116   default: assert(0);
00117   }
00118   return ret;
00119 }
00120 
00121 // ---------------------------------------------------------------------- 
00122 void Kernel3d_MPI::homogeneousDeg(vector<double>& degVec) const
00123 {
00124   switch(_kernelType) {
00125   case KNL_LAP_S_U: degVec.resize(1); degVec[0]=1; break;
00126   case KNL_STK_F_U: degVec.resize(4); degVec[0]=1; degVec[1]=1; degVec[2]=1; degVec[3]=2; break;
00127   case KNL_NAV_S_U: degVec.resize(3); degVec[0]=1; degVec[1]=1; degVec[2]=1; break;
00128   default: assert(0);
00129   }
00130   return;
00131 }
00132 
00133 // ----------------------------------------------------------------------
00134 /* This is the main function for Kernel3d_MPI.
00135  * Given source positions, normals and target positions, build a matrix inter
00136  * which will be used as a multiplier during the multiplication phase.
00137  * This functino is just a set of cases, based on the kernelType variable.
00138  * Mathematical details on these kernels are available in the papers
00139  */
00140 #undef __FUNCT__
00141 #define __FUNCT__ "Kernel3d_MPI::kernel"
00142 int Kernel3d_MPI::buildKnlIntCtx(const DblNumMat& srcPos, const DblNumMat& srcNor, const DblNumMat& trgPos, DblNumMat& inter)
00143 {
00144   //begin
00145   int srcDOF = this->srcDOF();
00146   int trgDOF = this->trgDOF();
00147   pA(srcPos.m()==dim() && srcNor.m()==dim() && trgPos.m()==dim());
00148   pA(srcPos.n()==srcNor.n());
00149   pA(srcPos.n()*srcDOF == inter.n());
00150   pA(trgPos.n()*trgDOF == inter.m());
00151   /* Single-Layer Laplacian Position/Velocity */
00152   if(       _kernelType==KNL_LAP_S_U) {
00153          //----------------------------------------------------------------------------------
00154          //---------------------------------
00155          double OOFP = 1.0/(4.0*M_PI);
00156          for(int i=0; i<trgPos.n(); i++) {
00157                 for(int j=0; j<srcPos.n(); j++) {
00158                   double x = trgPos(0,i) - srcPos(0,j); double y = trgPos(1,i) - srcPos(1,j); double z = trgPos(2,i) - srcPos(2,j);
00159                   double r2 = x*x + y*y + z*z;                                                    double r = sqrt(r2);
00160                   if(r<_mindif) {
00161                          inter(i,j) = 0;
00162                   } else {
00163                          inter(i,j) = OOFP / r;
00164                   }
00165                 }
00166          }
00167   }
00168   /* Double-Layer Laplacian Position/Velocity */
00169   else if(_kernelType==KNL_LAP_D_U) {
00170          //---------------------------------
00171          double OOFP = 1.0/(4.0*M_PI);
00172          for(int i=0; i<trgPos.n(); i++)
00173                 for(int j=0; j<srcPos.n(); j++) {
00174                   double x = trgPos(0,i) - srcPos(0,j); double y = trgPos(1,i) - srcPos(1,j); double z = trgPos(2,i) - srcPos(2,j);
00175                   double r2 = x*x + y*y + z*z;                    double r = sqrt(r2);
00176                   if(r<_mindif) {
00177                          for(int t=0;t<trgDOF;t++) for(int s=0;s<srcDOF;s++) { inter(i*trgDOF+t, j*srcDOF+s) = 0.0; }
00178                   } else {
00179                          double nx = srcNor(0,j);                                 double ny = srcNor(1,j);                                double nz = srcNor(2,j);
00180                          double rn = x*nx + y*ny + z*nz;
00181                          double r3 = r2*r;
00182                          inter(i,j) = - OOFP / r3 * rn;
00183                   }
00184                 }
00185   }
00186   /* Laplacian Identity Kernel */
00187   else if(_kernelType==KNL_LAP_I) {
00188          //---------------------------------
00189          for(int i=0; i<trgPos.n(); i++)
00190                 for(int j=0; j<srcPos.n(); j++)
00191                   inter(i,j) = 1;
00192   }
00193   /* Stokes - Fmm3d velocity kernel */
00194   else if(_kernelType==KNL_STK_F_U) {
00195          //----------------------------------------------------------------------------------
00196          //---------------------------------
00197          pA(_coefs.size()>=1);
00198          double mu = _coefs[0];
00199          double OOFP = 1.0/(4.0*M_PI);
00200          double OOEP = 1.0/(8.0*M_PI);
00201          double oomu = 1.0/mu;
00202          for(int i=0; i<trgPos.n(); i++)
00203                 for(int j=0; j<srcPos.n(); j++) {
00204                   double x = trgPos(0,i) - srcPos(0,j); double y = trgPos(1,i) - srcPos(1,j); double z = trgPos(2,i) - srcPos(2,j);
00205                   double r2 = x*x + y*y + z*z;                    double r = sqrt(r2);
00206                   if(r<_mindif) {
00207                          for(int t=0;t<trgDOF;t++) for(int s=0;s<srcDOF;s++) { inter(i*trgDOF+t, j*srcDOF+s) = 0.0; }
00208                   } else {
00209                          double r3 = r2*r;
00210                          double G = oomu * OOEP / r;
00211                          double H = oomu * OOEP / r3;
00212                          double A = OOFP / r3;
00213                          int is = i*3;                   int js = j*4;
00214                          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;
00215                          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;
00216                          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;
00217                   }
00218                 }
00219   }
00220   /* Stokes Single-Layer Position/Velocity kernel */
00221   else if(_kernelType==KNL_STK_S_U) {
00222          //---------------------------------
00223          pA(_coefs.size()>=1);
00224          double mu = _coefs[0];
00225          double OOEP = 1.0/(8.0*M_PI);
00226          double oomu = 1.0/mu;
00227          for(int i=0; i<trgPos.n(); i++)
00228                 for(int j=0; j<srcPos.n(); j++) {
00229                   double x = trgPos(0,i) - srcPos(0,j); double y = trgPos(1,i) - srcPos(1,j); double z = trgPos(2,i) - srcPos(2,j);
00230                   double r2 = x*x + y*y + z*z;                    double r = sqrt(r2);
00231                   if(r<_mindif) {
00232                          for(int t=0;t<trgDOF;t++) for(int s=0;s<srcDOF;s++) { inter(i*trgDOF+t, j*srcDOF+s) = 0.0; }
00233                   } else {
00234                          double r3 = r2*r;
00235                          double G = OOEP / r;
00236                          double H = OOEP / r3;
00237                          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);
00238                          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);
00239                          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);
00240                   }
00241                 }
00242   }
00243   /* Stokes Single-Layer Pressure kernel */
00244   else if(_kernelType==KNL_STK_S_P) {
00245          //---------------------------------
00246          pA(_coefs.size()>=1);
00247          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); double y = trgPos(1,i) - srcPos(1,j); double z = trgPos(2,i) - srcPos(2,j);
00252                   double r2 = x*x + y*y + z*z;                    double r = sqrt(r2);
00253                   if(r<_mindif) {
00254                          for(int t=0;t<trgDOF;t++) for(int s=0;s<srcDOF;s++) { inter(i*trgDOF+t, j*srcDOF+s) = 0.0; }
00255                   } else {
00256                          double r3 = r2*r;
00257                          inter(i  ,j*3  ) = OOFP*x/r3;                   inter(i  ,j*3+1) = OOFP*y/r3;                   inter(i  ,j*3+2) = OOFP*z/r3;
00258                   }
00259                 }
00260   }
00261   /* Stokes Double-Layer Position/Velocity kernel */
00262   else if(_kernelType==KNL_STK_D_U) {
00263          //---------------------------------
00264          pA(_coefs.size()>=1);
00265          double mu = _coefs[0];
00266          double SOEP = 6.0/(8.0*M_PI);
00267          for(int i=0; i<trgPos.n(); i++)
00268                 for(int j=0; j<srcPos.n(); j++) {
00269                   double x = trgPos(0,i) - srcPos(0,j); double y = trgPos(1,i) - srcPos(1,j); double z = trgPos(2,i) - srcPos(2,j);
00270                   double r2 = x*x + y*y + z*z;                    double r = sqrt(r2);
00271                   if(r<_mindif) {
00272                          for(int t=0;t<trgDOF;t++) for(int s=0;s<srcDOF;s++) { inter(i*trgDOF+t, j*srcDOF+s) = 0.0; }
00273                   } else {
00274                          double nx = srcNor(0,j);                                 double ny = srcNor(1,j);                                double nz = srcNor(2,j);
00275                          double rn = x*nx + y*ny + z*nz;
00276                          double r5 = r2*r2*r;
00277                          double C = - SOEP / r5;
00278                          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;
00279                          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;
00280                          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;
00281                   }
00282                 }
00283   }
00284   /* Stokes Double-Layer Pressure kernel */
00285   else if(_kernelType==KNL_STK_D_P) {
00286          //---------------------------------
00287          pA(_coefs.size()>=1);
00288          double mu = _coefs[0];
00289          double OOTP = 1.0/(2.0*M_PI);
00290          double coef = mu*OOTP;
00291          for(int i=0; i<trgPos.n(); i++)
00292                 for(int j=0; j<srcPos.n(); j++) {
00293                   double x = trgPos(0,i) - srcPos(0,j); double y = trgPos(1,i) - srcPos(1,j); double z = trgPos(2,i) - srcPos(2,j);
00294                   double r2 = x*x + y*y + z*z;                    double r = sqrt(r2);
00295                   if(r<_mindif) {
00296                          for(int t=0;t<trgDOF;t++) for(int s=0;s<srcDOF;s++) { inter(i*trgDOF+t, j*srcDOF+s) = 0.0; }
00297                   } else {
00298                          double nx = srcNor(0,j);                                 double ny = srcNor(1,j);                                double nz = srcNor(2,j);
00299                          double rn = x*nx + y*ny + z*nz;
00300                          double r3 = r2*r;
00301                          double r5 = r3*r2;
00302                          int is = i;                     int js = j*3;
00303                          inter(is  ,js  ) = coef*(nx/r3 - 3*rn*x/r5);
00304                          inter(is  ,js+1) = coef*(ny/r3 - 3*rn*y/r5);
00305                          inter(is  ,js+2) = coef*(nz/r3 - 3*rn*z/r5);
00306                   }
00307                 }
00308   }
00309   else if(_kernelType==KNL_STK_R_U) {
00310          //---------------------------------
00311          pA(_coefs.size()>=1);
00312          double mu = _coefs[0];
00313          for(int i=0; i<trgPos.n(); i++)
00314                 for(int j=0; j<srcPos.n(); j++) {
00315                   double x = trgPos(0,i) - srcPos(0,j); double y = trgPos(1,i) - srcPos(1,j); double z = trgPos(2,i) - srcPos(2,j);
00316                   double r2 = x*x + y*y + z*z;                    double r = sqrt(r2);
00317                   if(r<_mindif) {
00318                          for(int t=0;t<trgDOF;t++) for(int s=0;s<srcDOF;s++) { inter(i*trgDOF+t, j*srcDOF+s) = 0.0; }
00319                   } else {
00320                          double r3 = r2*r;
00321                          double coef = 1.0/(8.0*M_PI*mu)/r3;
00322                          inter(i*3,   j*3)   = coef*0;                   inter(i*3,   j*3+1) = coef*z;                   inter(i*3,   j*3+2) = coef*(-y);
00323                          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;
00324                          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;
00325                   }
00326                 }
00327   } else if(_kernelType==KNL_STK_R_P) {
00328          //---------------------------------
00329          pA(_coefs.size()>=1);
00330          double mu = _coefs[0];
00331          for(int i=0; i<trgPos.n(); i++)
00332                 for(int j=0; j<srcPos.n(); j++) {
00333                   inter(i,j*3  ) = 0;
00334                   inter(i,j*3+1) = 0;
00335                   inter(i,j*3+2) = 0;
00336                 }
00337   } else if(_kernelType==KNL_STK_I) {
00338          //---------------------------------
00339          pA(_coefs.size()>=1);
00340          double mu = _coefs[0];
00341          for(int i=0; i<trgPos.n(); i++)
00342                 for(int j=0; j<srcPos.n(); j++) {
00343                   inter(i*3,   j*3  ) = 1;               inter(i*3,   j*3+1) = 0;                inter(i*3,   j*3+2) = 0;
00344                   inter(i*3+1, j*3  ) = 0;               inter(i*3+1, j*3+1) = 1;                inter(i*3+1, j*3+2) = 0;
00345                   inter(i*3+2, j*3  ) = 0;               inter(i*3+2, j*3+1) = 0;                inter(i*3+2, j*3+2) = 1;
00346                 }
00347   }
00348   /* levi-civita tensor */
00349   else if(_kernelType==KNL_STK_E) {
00350          //---------------------------------
00351          pA(_coefs.size()>=1);
00352          double mu = _coefs[0];
00353          for(int i=0; i<trgPos.n(); i++)
00354                 for(int j=0; j<srcPos.n(); j++) {
00355                   double x = trgPos(0,i) - srcPos(0,j); double y = trgPos(1,i) - srcPos(1,j); double z = trgPos(2,i) - srcPos(2,j);
00356                   inter(i*3,   j*3  ) = 0;               inter(i*3,   j*3+1) = z;                inter(i*3,   j*3+2) = -y;
00357                   inter(i*3+1, j*3  ) = -z;      inter(i*3+1, j*3+1) = 0;                inter(i*3+1, j*3+2) = x;
00358                   inter(i*3+2, j*3  ) = y;               inter(i*3+2, j*3+1) = -x;               inter(i*3+2, j*3+2) = 0;
00359                 }
00360   }
00361   /* Navier-Stokes Single-Layer Position/Velocity Kernel */
00362   else if(_kernelType==KNL_NAV_S_U) {
00363          //----------------------------------------------------------------------------------
00364          //---------------------------------
00365          pA(_coefs.size()>=2);
00366          double mu = _coefs[0];  double ve = _coefs[1];
00367          double sc1 = (3.0-4.0*ve)/(16.0*M_PI*(1.0-ve));
00368          double sc2 = 1.0/(16.0*M_PI*(1.0-ve));
00369          double oomu = 1.0/mu;
00370          for(int i=0; i<trgPos.n(); i++)
00371                 for(int j=0; j<srcPos.n(); j++) {
00372                   double x = trgPos(0,i) - srcPos(0,j); double y = trgPos(1,i) - srcPos(1,j); double z = trgPos(2,i) - srcPos(2,j);
00373                   double r2 = x*x + y*y + z*z;                    double r = sqrt(r2);
00374                   if(r<_mindif) {
00375                          for(int t=0;t<trgDOF;t++) for(int s=0;s<srcDOF;s++) { inter(i*trgDOF+t, j*srcDOF+s) = 0.0; }
00376                   } else {
00377                          double r3 = r2*r;
00378                          double G = sc1 / r;
00379                          double H = sc2 / r3;
00380                          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);
00381                          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);
00382                          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);
00383                   }
00384                 }
00385   }
00386   /* Navier-Stokes Double-Layer Position/Velocity Kernel */
00387   else if(_kernelType==KNL_NAV_D_U) {
00388          //---------------------------------
00389          pA(_coefs.size()>=2);
00390          double mu = _coefs[0];  double ve = _coefs[1];
00391          double dc1 = -(1-2.0*ve)/(8.0*M_PI*(1.0-ve));
00392          double dc2 =  (1-2.0*ve)/(8.0*M_PI*(1.0-ve));
00393          double dc3 = -3.0/(8.0*M_PI*(1.0-ve));
00394          for(int i=0; i<trgPos.n(); i++)
00395                 for(int j=0; j<srcPos.n(); j++) {
00396                   double x = trgPos(0,i) - srcPos(0,j); double y = trgPos(1,i) - srcPos(1,j); double z = trgPos(2,i) - srcPos(2,j);
00397                   double r2 = x*x + y*y + z*z;                    double r = sqrt(r2);
00398                   if(r<_mindif) {
00399                          for(int t=0;t<trgDOF;t++) for(int s=0;s<srcDOF;s++) { inter(i*trgDOF+t, j*srcDOF+s) = 0.0; }
00400                   } else {
00401                          double nx = srcNor(0,j);                                 double ny = srcNor(1,j);                                double nz = srcNor(2,j);
00402                          double rn = x*nx + y*ny + z*nz;
00403                          double r3 = r2*r;
00404                          double r5 = r3*r2;
00405                          double A = dc1 / r3;                    double B = dc2 / r3;                    double C = dc3 / r5;
00406                          double rx = x;                  double ry = y;                  double rz = z;
00407                          //&&&&&&&
00408                          inter(i*3,   j*3)   = A*(rn+nx*rx) + B*(rx*nx) + C*rn*rx*rx;
00409                          inter(i*3,   j*3+1) = A*(   nx*ry) + B*(rx*ny) + C*rn*rx*ry;
00410                          inter(i*3,   j*3+2) = A*(   nx*rz) + B*(rx*nz) + C*rn*rx*rz;
00411                          inter(i*3+1, j*3)   = A*(   ny*rx) + B*(ry*nx) + C*rn*ry*rx;
00412                          inter(i*3+1, j*3+1) = A*(rn+ny*ry) + B*(ry*ny) + C*rn*ry*ry;
00413                          inter(i*3+1, j*3+2) = A*(   ny*rz) + B*(ry*nz) + C*rn*ry*rz;
00414                          inter(i*3+2, j*3)   = A*(   nz*rx) + B*(rz*nx) + C*rn*rz*rx;
00415                          inter(i*3+2, j*3+1) = A*(   nz*ry) + B*(rz*ny) + C*rn*rz*ry;
00416                          inter(i*3+2, j*3+2) = A*(rn+nz*rz) + B*(rz*nz) + C*rn*rz*rz;
00417                   }
00418                 }
00419   } else if(_kernelType==KNL_NAV_R_U) {
00420          //---------------------------------
00421          pA(_coefs.size()>=2);
00422          double mu = _coefs[0];  double ve = _coefs[1];
00423          for(int i=0; i<trgPos.n(); i++)
00424                 for(int j=0; j<srcPos.n(); j++) {
00425                   double x = trgPos(0,i) - srcPos(0,j); double y = trgPos(1,i) - srcPos(1,j); double z = trgPos(2,i) - srcPos(2,j);
00426                   double r2 = x*x + y*y + z*z;                    double r = sqrt(r2);
00427                   if(r<_mindif) {
00428                          for(int t=0;t<trgDOF;t++) for(int s=0;s<srcDOF;s++) { inter(i*trgDOF+t, j*srcDOF+s) = 0.0; }
00429                   } else {
00430                          double r3 = r2*r;
00431                          double coef = 1.0/(8.0*M_PI*mu)/r3;
00432                          inter(i*3,   j*3)   = coef*0;                   inter(i*3,   j*3+1) = coef*z;                   inter(i*3,   j*3+2) = coef*(-y);
00433                          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;
00434                          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;
00435                   }
00436                 }
00437   } else if(_kernelType==KNL_NAV_I) {
00438          //---------------------------------
00439          pA(_coefs.size()>=2);
00440          double mu = _coefs[0];  double ve = _coefs[1];
00441          for(int i=0; i<trgPos.n(); i++)
00442                 for(int j=0; j<srcPos.n(); j++) {
00443                   inter(i*3,   j*3  ) = 1;               inter(i*3,   j*3+1) = 0;                inter(i*3,   j*3+2) = 0;
00444                   inter(i*3+1, j*3  ) = 0;               inter(i*3+1, j*3+1) = 1;                inter(i*3+1, j*3+2) = 0;
00445                   inter(i*3+2, j*3  ) = 0;               inter(i*3+2, j*3+1) = 0;                inter(i*3+2, j*3+2) = 1;
00446                 }
00447   } else if(_kernelType==KNL_NAV_E) {
00448          //---------------------------------
00449          pA(_coefs.size()>=2);
00450          double mu = _coefs[0];  double ve = _coefs[1];
00451          for(int i=0; i<trgPos.n(); i++)
00452                 for(int j=0; j<srcPos.n(); j++) {
00453                   double x = trgPos(0,i) - srcPos(0,j); double y = trgPos(1,i) - srcPos(1,j); double z = trgPos(2,i) - srcPos(2,j);
00454                   inter(i*3,   j*3  ) = 0;               inter(i*3,   j*3+1) = z;                inter(i*3,   j*3+2) = -y;
00455                   inter(i*3+1, j*3  ) = -z;      inter(i*3+1, j*3+1) = 0;                inter(i*3+1, j*3+2) = x;
00456                   inter(i*3+2, j*3  ) = y;               inter(i*3+2, j*3+1) = -x;               inter(i*3+2, j*3+2) = 0;
00457                 }
00458   } else if(_kernelType==KNL_ERR) {
00459          //---------------------------------
00460          pA(0);
00461   }
00462   return(0);
00463 }
00464 
00465 
00466 
00467 

Generated on Sun Dec 4 21:12:40 2005 for fmm3d_mpi by  doxygen 1.4.5