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<Triangle::prm> & prmTri ,
               const int ScalfmmNbThreads )
        : prmTetra(prmTet), prmTriangle(prmTri),
          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 < (FSize)msh.magNode.size(); ++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<Triangle::Tri, Triangle::NPI>(msh.tri, msh.magTri, idxPart, msh.c);
 
        srcDen.resize( msh.magTri.size()*Triangle::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);
        if (!FIRST_ORDER)
            { demag(Nodes::get_v<Nodes::NEXT>, Nodes::set_phiv, msh); }
        }
 
    std::vector<double> corr;
 
    const std::vector<Tetra::prm> &prmTetra;
 
    const std::vector<Triangle::prm> &prmTriangle;
 
private:
    std::vector<double> srcDen;
 
    OctreeClass tree;
 
    KernelClass kernels;
 
    double norm;
 
    template<class T, const int NPI>
    void insertCharges(const std::vector<T> &container, const std::vector<int> &idxContainer,
                       FSize &idx, const Eigen::Ref<const Eigen::Vector3d> c)
        {
        std::for_each(idxContainer.begin(), idxContainer.end(),
                      [this,&container, c, &idx](const int idxElem)
                      {
                      const T &elem = container[idxElem];
                      Eigen::Matrix<double,Nodes::DIM,NPI> gauss = elem.getPtGauss();
 
                      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(const std::function<const Eigen::Vector3d(const 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];
                Eigen::Matrix<double,Tetra::NPI,1> result =
                        t.charges(prmTetra[t.idxPrm].Ms, getter);
                for(int i=0;i<Tetra::NPI;i++)
                    { srcDen[nsrc+i] = result(i); }
                nsrc += Tetra::NPI;
                });
        std::fill(corr.begin(),corr.end(),0);
        std::for_each(msh.magTri.begin(),msh.magTri.end(),[this, &msh, &getter, &nsrc](const int idx)
                {
                Triangle::Tri &f = msh.tri[idx];
                Eigen::Matrix<double,Triangle::NPI,1> result = f.charges(f.dMs, getter);
                for(int i=0;i<Triangle::NPI;i++)
                    { srcDen[nsrc+i] = result(i); }
                nsrc += Triangle::NPI;
                f.correctionCharges(getter,result,corr);
                });
        }
 
    void demag(const std::function<const Eigen::Vector3d(const Nodes::Node&)>& getter,
               const 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