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), NOD(msh.getNbNodes()),
          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 < NOD; ++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, idxPart, msh.c);
        insertCharges<Facette::Fac, Facette::NPI>(msh.fac, idxPart, msh.c);
 
        srcDen.resize( msh.getNbFacs()*Facette::NPI + msh.getNbTets()*Tetra::NPI );
        corr.resize(NOD);
        }
 
    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;
 
    std::vector<Tetra::prm> prmTetra;
 
    std::vector<Facette::prm> prmFacette;
 
private:
    const int NOD; 
    OctreeClass tree;    
    KernelClass kernels; 
    double norm; 
    template<class T, const int NPI>
    void insertCharges(std::vector<T> const &container, FSize &idx, Eigen::Ref<Eigen::Vector3d> const c)
        {
        std::for_each(container.begin(), container.end(),
                      [this, c, &idx](T const &elem)
                      {
                          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);
                              }
                      });  // end for_each
        }
 
    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.tet.begin(), msh.tet.end(),
                      [this, getter, &nsrc](Tetra::Tet const &tet)
                          { tet.charges(prmTetra[tet.idxPrm], getter, srcDen, nsrc); });
        std::fill(corr.begin(),corr.end(),0);
        std::for_each(msh.fac.begin(), msh.fac.end(),
                      [this, getter, &nsrc](Facette::Fac const &fac)
                          { fac.charges(prmFacette[fac.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);
 
        // reset potentials and forces - physicalValues[idxPart] = Q
        tree.forEachLeaf(
                [this](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] - NOD];
                        }
 
                    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