FP2PTensorialKij.hpp 15.4 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
// ===================================================================================
#ifndef FP2P_TENSORIALKIJ_HPP
#define FP2P_TENSORIALKIJ_HPP

namespace FP2P {

//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
// Tensorial Matrix Kernels: K_IJ / p_i=\sum_j K_{ij} w_j
//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////

/**
   * @brief MutualParticlesKIJ
   * @param sourceX
   * @param sourceY
   * @param sourceZ
   * @param sourcePhysicalValue
   * @param sourceForceX
   * @param sourceForceY
   * @param sourceForceZ
   * @param sourcePotential
   * @param targetX
   * @param targetY
   * @param targetZ
   * @param targetPhysicalValue
   * @param targetForceX
   * @param targetForceY
   * @param targetForceZ
   * @param targetPotential
   */
template<class FReal, typename MatrixKernelClass>
inline void MutualParticlesKIJ(const FReal targetX,const FReal targetY,const FReal targetZ, const FReal* targetPhysicalValue,
                               FReal* targetForceX, FReal* targetForceY, FReal* targetForceZ, FReal* targetPotential,
                               const FReal sourceX,const FReal sourceY,const FReal sourceZ, const FReal* sourcePhysicalValue,
                               FReal* sourceForceX, FReal* sourceForceY, FReal* sourceForceZ, FReal* sourcePotential,
                               const MatrixKernelClass *const MatrixKernel){

    // get information on tensorial aspect of matrix kernel
    const int ncmp = MatrixKernelClass::NCMP;    
    const int npv = MatrixKernelClass::NPV;    
    const int npot = MatrixKernelClass::NPOT;    

    // evaluate kernel and its partial derivatives
    const FPoint<FReal> sourcePoint(sourceX,sourceY,sourceZ);
    const FPoint<FReal> targetPoint(targetX,targetY,targetZ);
    FReal Kxy[ncmp];
    FReal dKxy[ncmp][3];
    MatrixKernel->evaluateBlockAndDerivative(targetPoint,sourcePoint,Kxy,dKxy);
    const FReal mutual_coeff = MatrixKernel->getMutualCoefficient(); // 1 if symmetric; -1 if antisymmetric

    for(unsigned int i = 0 ; i < npot ; ++i){
        for(unsigned int j = 0 ; j < npv ; ++j){

            // update component to be applied
            const int d = MatrixKernel->getPosition(i*npv+j);

            // forces prefactor
            const FReal coef = (targetPhysicalValue[j] * sourcePhysicalValue[j]);

            targetForceX[i] += dKxy[d][0] * coef;
            targetForceY[i] += dKxy[d][1] * coef;
            targetForceZ[i] += dKxy[d][2] * coef;
            targetPotential[i] += ( Kxy[d] * sourcePhysicalValue[j] );

            sourceForceX[i] -= dKxy[d][0] * coef;
            sourceForceY[i] -= dKxy[d][1] * coef;
            sourceForceZ[i] -= dKxy[d][2] * coef;
            sourcePotential[i] += ( mutual_coeff * Kxy[d] * targetPhysicalValue[j] );

        }// j
    }// i

}

//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
// Tensorial Matrix Kernels: K_IJ
// TODO: Implement SSE and AVX variants then move following FullMutualKIJ and FullRemoteKIJ to FP2P.hpp
//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////

/**
 * @brief FullMutualKIJ
 */
template <class FReal, class ContainerClass, typename MatrixKernelClass>
inline void FullMutualKIJ(ContainerClass* const FRestrict inTargets, ContainerClass* const inNeighbors[],
                          const int limiteNeighbors, const MatrixKernelClass *const MatrixKernel){

    // get information on tensorial aspect of matrix kernel
    const int ncmp = MatrixKernelClass::NCMP;
    const int npv = MatrixKernelClass::NPV;    
    const int npot = MatrixKernelClass::NPOT;

    // get info on targets
    const FSize nbParticlesTargets = inTargets->getNbParticles();
    const FReal*const targetsX = inTargets->getPositions()[0];
    const FReal*const targetsY = inTargets->getPositions()[1];
    const FReal*const targetsZ = inTargets->getPositions()[2];

    for(FSize idxNeighbors = 0 ; idxNeighbors < limiteNeighbors ; ++idxNeighbors){
        for(FSize idxTarget = 0 ; idxTarget < nbParticlesTargets ; ++idxTarget){
            if( inNeighbors[idxNeighbors] ){
                const FSize nbParticlesSources = inNeighbors[idxNeighbors]->getNbParticles();
                const FReal*const sourcesX = inNeighbors[idxNeighbors]->getPositions()[0];
                const FReal*const sourcesY = inNeighbors[idxNeighbors]->getPositions()[1];
                const FReal*const sourcesZ = inNeighbors[idxNeighbors]->getPositions()[2];

                for(FSize idxSource = 0 ; idxSource < nbParticlesSources ; ++idxSource){

                    // evaluate kernel and its partial derivatives
                    const FPoint<FReal> sourcePoint(sourcesX[idxSource],sourcesY[idxSource],sourcesZ[idxSource]);
                    const FPoint<FReal> targetPoint(targetsX[idxTarget],targetsY[idxTarget],targetsZ[idxTarget]);
                    FReal Kxy[ncmp];
                    FReal dKxy[ncmp][3];
                    MatrixKernel->evaluateBlockAndDerivative(targetPoint,sourcePoint,Kxy,dKxy);
                    const FReal mutual_coeff = MatrixKernel->getMutualCoefficient(); // 1 if symmetric; -1 if antisymmetric

                    for(unsigned int i = 0 ; i < npot ; ++i){
                        FReal*const targetsPotentials = inTargets->getPotentials(i);
                        FReal*const targetsForcesX = inTargets->getForcesX(i);
                        FReal*const targetsForcesY = inTargets->getForcesY(i);
                        FReal*const targetsForcesZ = inTargets->getForcesZ(i);
                        FReal*const sourcesPotentials = inNeighbors[idxNeighbors]->getPotentials(i);
                        FReal*const sourcesForcesX = inNeighbors[idxNeighbors]->getForcesX(i);
                        FReal*const sourcesForcesY = inNeighbors[idxNeighbors]->getForcesY(i);
                        FReal*const sourcesForcesZ = inNeighbors[idxNeighbors]->getForcesZ(i);

                        for(unsigned int j = 0 ; j < npv ; ++j){
                            const FReal*const targetsPhysicalValues = inTargets->getPhysicalValues(j);
                            const FReal*const sourcesPhysicalValues = inNeighbors[idxNeighbors]->getPhysicalValues(j);

                            // update component to be applied
                            const int d = MatrixKernel->getPosition(i*npv+j);

                            // forces prefactor
                            FReal coef = (targetsPhysicalValues[idxTarget] * sourcesPhysicalValues[idxSource]);

                            targetsForcesX[idxTarget] += dKxy[d][0] * coef;
                            targetsForcesY[idxTarget] += dKxy[d][1] * coef;
                            targetsForcesZ[idxTarget] += dKxy[d][2] * coef;
                            targetsPotentials[idxTarget] += ( Kxy[d] * sourcesPhysicalValues[idxSource] );

                            sourcesForcesX[idxSource] -= dKxy[d][0] * coef;
                            sourcesForcesY[idxSource] -= dKxy[d][1] * coef;
                            sourcesForcesZ[idxSource] -= dKxy[d][2] * coef;
                            sourcesPotentials[idxSource] += mutual_coeff * Kxy[d] * targetsPhysicalValues[idxTarget];

                        }// j
                    }// i
                }
            }
        }
    }
}

template <class FReal, class ContainerClass, typename MatrixKernelClass>
inline void InnerKIJ(ContainerClass* const FRestrict inTargets, const MatrixKernelClass *const MatrixKernel){

    // get information on tensorial aspect of matrix kernel
    const int ncmp = MatrixKernelClass::NCMP;
    const int npv = MatrixKernelClass::NPV;    
    const int npot = MatrixKernelClass::NPOT;

    // get info on targets
    const FSize nbParticlesTargets = inTargets->getNbParticles();
    const FReal*const targetsX = inTargets->getPositions()[0];
    const FReal*const targetsY = inTargets->getPositions()[1];
    const FReal*const targetsZ = inTargets->getPositions()[2];

    for(FSize idxTarget = 0 ; idxTarget < nbParticlesTargets ; ++idxTarget){
        for(FSize idxSource = idxTarget + 1 ; idxSource < nbParticlesTargets ; ++idxSource){

            // evaluate kernel and its partial derivatives
            const FPoint<FReal> sourcePoint(targetsX[idxSource],targetsY[idxSource],targetsZ[idxSource]);
            const FPoint<FReal> targetPoint(targetsX[idxTarget],targetsY[idxTarget],targetsZ[idxTarget]);
            FReal Kxy[ncmp];
            FReal dKxy[ncmp][3];
            MatrixKernel->evaluateBlockAndDerivative(targetPoint,sourcePoint,Kxy,dKxy);
            const FReal mutual_coeff = MatrixKernel->getMutualCoefficient(); // 1 if symmetric; -1 if antisymmetric

            for(unsigned int i = 0 ; i < npot ; ++i){
                FReal*const targetsPotentials = inTargets->getPotentials(i);
                FReal*const targetsForcesX = inTargets->getForcesX(i);
                FReal*const targetsForcesY = inTargets->getForcesY(i);
                FReal*const targetsForcesZ = inTargets->getForcesZ(i);

                for(unsigned int j = 0 ; j < npv ; ++j){
                    const FReal*const targetsPhysicalValues = inTargets->getPhysicalValues(j);

                    // update component to be applied
                    const int d = MatrixKernel->getPosition(i*npv+j);

                    // forces prefactor
                    const FReal coef = (targetsPhysicalValues[idxTarget] * targetsPhysicalValues[idxSource]);

                    targetsForcesX[idxTarget] += dKxy[d][0] * coef;
                    targetsForcesY[idxTarget] += dKxy[d][1] * coef;
                    targetsForcesZ[idxTarget] += dKxy[d][2] * coef;
                    targetsPotentials[idxTarget] += ( Kxy[d] * targetsPhysicalValues[idxSource] );

                    targetsForcesX[idxSource] -= dKxy[d][0] * coef;
                    targetsForcesY[idxSource] -= dKxy[d][1] * coef;
                    targetsForcesZ[idxSource] -= dKxy[d][2] * coef;
                    targetsPotentials[idxSource] += mutual_coeff * Kxy[d] * targetsPhysicalValues[idxTarget];
                }// j
            }// i

        }
    }
}

/**
   * @brief FullRemoteKIJ
   */
template <class FReal, class ContainerClass, typename MatrixKernelClass>
inline void FullRemoteKIJ(ContainerClass* const FRestrict inTargets, const ContainerClass* const inNeighbors[],
                          const int limiteNeighbors, const MatrixKernelClass *const MatrixKernel){

    // get information on tensorial aspect of matrix kernel
    const int ncmp = MatrixKernelClass::NCMP;
    const int npv = MatrixKernelClass::NPV;    
    const int npot = MatrixKernelClass::NPOT;

    // get info on targets
    const FSize nbParticlesTargets = inTargets->getNbParticles();
    const FReal*const targetsX = inTargets->getPositions()[0];
    const FReal*const targetsY = inTargets->getPositions()[1];
    const FReal*const targetsZ = inTargets->getPositions()[2];

    for(FSize idxNeighbors = 0 ; idxNeighbors < limiteNeighbors ; ++idxNeighbors){
        for(FSize idxTarget = 0 ; idxTarget < nbParticlesTargets ; ++idxTarget){
            if( inNeighbors[idxNeighbors] ){
                const FSize nbParticlesSources = inNeighbors[idxNeighbors]->getNbParticles();
                const FReal*const sourcesX = inNeighbors[idxNeighbors]->getPositions()[0];
                const FReal*const sourcesY = inNeighbors[idxNeighbors]->getPositions()[1];
                const FReal*const sourcesZ = inNeighbors[idxNeighbors]->getPositions()[2];

                for(FSize idxSource = 0 ; idxSource < nbParticlesSources ; ++idxSource){

                    // evaluate kernel and its partial derivatives
                    const FPoint<FReal> sourcePoint(sourcesX[idxSource],sourcesY[idxSource],sourcesZ[idxSource]);
                    const FPoint<FReal> targetPoint(targetsX[idxTarget],targetsY[idxTarget],targetsZ[idxTarget]);
                    FReal Kxy[ncmp];
                    FReal dKxy[ncmp][3];
                    MatrixKernel->evaluateBlockAndDerivative(targetPoint,sourcePoint,Kxy,dKxy);

                    for(unsigned int i = 0 ; i < npot ; ++i){
                        FReal*const targetsPotentials = inTargets->getPotentials(i);
                        FReal*const targetsForcesX = inTargets->getForcesX(i);
                        FReal*const targetsForcesY = inTargets->getForcesY(i);
                        FReal*const targetsForcesZ = inTargets->getForcesZ(i);

                        for(unsigned int j = 0 ; j < npv ; ++j){
                            const FReal*const targetsPhysicalValues = inTargets->getPhysicalValues(j);
                            const FReal*const sourcesPhysicalValues = inNeighbors[idxNeighbors]->getPhysicalValues(j);

                            // update component to be applied
                            const int d = MatrixKernel->getPosition(i*npv+j);

                            // forces prefactor
                            const FReal coef = (targetsPhysicalValues[idxTarget] * sourcesPhysicalValues[idxSource]);

                            targetsForcesX[idxTarget] += dKxy[d][0] * coef;
                            targetsForcesY[idxTarget] += dKxy[d][1] * coef;
                            targetsForcesZ[idxTarget] += dKxy[d][2] * coef;
                            targetsPotentials[idxTarget] += ( Kxy[d] * sourcesPhysicalValues[idxSource] );

                        }// j
                    }// i

                }
            }
        }
    }
}



} // End namespace


#endif // FP2P_TENSORIALKIJ_HPP