fmm_demag.h

Go to the documentation of this file.

#ifndef FMM_DEMAG_H
#define FMM_DEMAG_H
 
#include "Components/FParticleType.hpp"
#include "Components/FTypedLeaf.hpp"
#include "Containers/FOctree.hpp"
#include "Core/FFmmAlgorithmThreadTsm.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
#include "Kernels/Rotation/FRotationCell.hpp"
#include "Kernels/Rotation/FRotationKernel.hpp"
 
#include "mesh.h"
 
namespace scal_fmm
    {
const int P = 9;              
const int NbLevels = 6;       
const int SizeSubLevels = 3;  
typedef double FReal; 
typedef FTypedRotationCell<FReal, P>
        CellClass; 
typedef FP2PParticleContainerIndexed<FReal>
        ContainerClass; 
typedef FTypedLeaf<FReal, ContainerClass>
        LeafClass; 
typedef FOctree<FReal, CellClass, ContainerClass, LeafClass>
        OctreeClass; 
typedef FRotationKernel<FReal, CellClass, ContainerClass, P>
        KernelClass; 
typedef FFmmAlgorithmThreadTsm<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass>
        FmmClass; 
const double boxWidth = 2.01;              
const FPoint<FReal> boxCenter(0., 0., 0.); 
class fmm
    {
public:
    inline fmm(Mesh::mesh &msh ,
               std::vector<Tetra::prm> & prmTet ,
               std::vector<Facette::prm> & prmFac ,
               const int ScalfmmNbThreads )
        : prmTetra(prmTet), prmFacette(prmFac),
          tree(NbLevels, SizeSubLevels, boxWidth, boxCenter), kernels(NbLevels, boxWidth, boxCenter)
        {
        omp_set_num_threads(ScalfmmNbThreads);
        norm = 2. / msh.l.maxCoeff();
 
        FSize idxPart = 0;
        for (idxPart = 0; idxPart < msh.getNbNodes(); ++idxPart)
            {
            if (msh.magNode[idxPart])
                {
                Eigen::Vector3d pTarget = norm*(msh.getNode_p(idxPart) - msh.c);
                tree.insert(FPoint<FReal>(pTarget.x(), pTarget.y(), pTarget.z()),
                            FParticleType::FParticleTypeTarget, idxPart);
                }
            }
 
        insertCharges<Tetra::Tet, Tetra::NPI>(msh.tet, msh.magTet, idxPart, msh.c);
        insertCharges<Facette::Fac, Facette::NPI>(msh.fac, msh.magFac, idxPart, msh.c);
 
        srcDen.resize( msh.magFac.size()*Facette::NPI + msh.magTet.size()*Tetra::NPI );
        corr.resize(msh.magNode.size());
        }
 
    void calc_demag(Mesh::mesh &msh )
        {
        demag(Nodes::get_u<Nodes::NEXT>, Nodes::set_phi, msh);
        demag(Nodes::get_v<Nodes::NEXT>, Nodes::set_phiv, msh);
        }
 
    std::vector<double> srcDen;
    
    std::vector<double> corr;
 
    const std::vector<Tetra::prm> &prmTetra;
 
    const std::vector<Facette::prm> &prmFacette;
 
private:
    OctreeClass tree;
 
    KernelClass kernels;
 
    double norm;
 
    template<class T, const int NPI>
    void insertCharges(std::vector<T> const &container, std::vector<int> const &idxContainer,
                       FSize &idx, Eigen::Ref<Eigen::Vector3d> const c)
        {
        std::for_each(idxContainer.begin(), idxContainer.end(),
                      [this,&container, c, &idx](int const &idxElem)
                      {
                      T const &elem = container[idxElem];
                      Eigen::Matrix<double,Nodes::DIM,NPI> gauss;
                      elem.getPtGauss(gauss);
 
                      for (int j = 0; j < NPI; j++, idx++)
                          {
                          double x = norm*(gauss(0,j) - c.x());
                          double y = norm*(gauss(1,j) - c.y());
                          double z = norm*(gauss(2,j) - c.z());
                          tree.insert(FPoint<FReal>(x,y,z), FParticleType::FParticleTypeSource, idx, 0.0);
                          }
                      });
        }
 
    void calc_charges(std::function<const Eigen::Vector3d(Nodes::Node)> getter, Mesh::mesh &msh)
        {
        int nsrc(0);
        std::fill(srcDen.begin(),srcDen.end(),0);
        std::for_each(msh.magTet.begin(),msh.magTet.end(),[this, &msh, getter, &nsrc](const int idx)
                {
                Tetra::Tet &t = msh.tet[idx];
                t.charges(prmTetra[t.idxPrm], getter, srcDen, nsrc);
                });
        std::fill(corr.begin(),corr.end(),0);
        std::for_each(msh.magFac.begin(),msh.magFac.end(),[this, &msh, getter, &nsrc](const int idx)
                {
                Facette::Fac &f = msh.fac[idx];
                f.charges(prmFacette[f.idxPrm], getter, srcDen, nsrc, corr);
                });
        }
 
    void demag(std::function<const Eigen::Vector3d(Nodes::Node)> getter,
               std::function<void(Nodes::Node &, const double)> setter, Mesh::mesh &msh)
        {
        FmmClass algo(&tree, &kernels);
        calc_charges(getter, msh);
 
        auto nbMagNodes = msh.magNode.size();
        // reset potentials and forces - physicalValues[idxPart] = Q
        tree.forEachLeaf(
                [this,&nbMagNodes](LeafClass *leaf)
                {
                    const int nbParticlesInLeaf = leaf->getSrc()->getNbParticles();
                    const auto &indexes = leaf->getSrc()->getIndexes();
                    FReal *const physicalValues = leaf->getSrc()->getPhysicalValues();
                    for (int idxPart = 0; idxPart < nbParticlesInLeaf; ++idxPart)
                        {
                        physicalValues[idxPart] = srcDen[indexes[idxPart] - nbMagNodes];
                        }
 
                    std::fill_n(leaf->getTargets()->getPotentials(),
                                leaf->getTargets()->getNbParticles(), 0);
                });
 
        tree.forEachCell([](CellClass *cell) { cell->resetToInitialState(); });
 
        algo.execute();
 
        tree.forEachLeaf(
                [this, &msh, setter](LeafClass *leaf)
                {
                    const FReal *const potentials = leaf->getTargets()->getPotentials();
                    const int nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
                    const auto &indexes = leaf->getTargets()->getIndexes();
                    for (int idxPart = 0; idxPart < nbParticlesInLeaf; ++idxPart)
                        {
                        const int indexPartOrig = indexes[idxPart];
                        msh.set(indexPartOrig, setter,
                                (potentials[idxPart] * norm + corr[indexPartOrig]) / (4 * M_PI));
                        }
                });
        }
    };  // end class fmm
 
    }  // namespace scal_fmm
#endif