fem.h

Go to the documentation of this file.

#ifndef fem_h
#define fem_h
 
#include <cassert>
#include <cmath>
#include <cstdlib>
#include <fstream>
#include <functional>
#include <iomanip>
#include <iostream>
#include <list>
#include <map>
#include <sstream>
#include <utility>
#include <vector>
 
#include <sys/times.h>
#include <time.h>
#include <unistd.h>
 
#include "mesh.h"
#include "fmm_demag.h"
#include "spinAccumulationSolver.h"
#include "linear_algebra.h"
#include "chronometer.h"
#include "ANN.h"
#include "feellgoodSettings.h"
 
const int NB_ENERGY_TERMS = 4;
 
enum ENERGY_TYPE
    {
    EXCHANGE = 0,
    ANISOTROPY = 1,
    DEMAG = 2,
    ZEEMAN = 3
    };
 
class Fem
    {
public:
    inline Fem(Settings &mySets, timing &t_prm) : msh(mySets)
        {
        vmax = 0.0;
        std::fill(E.begin(),E.end(),0);
        Etot0 = INFINITY;  // avoid "WARNING: energy increased" on first time step
        Etot = 0.0;
 
        recenter_mem = false;
        if (mySets.recenter)
            {
            if (mySets.verbose)
                {
                std::cout << "Approximate nearest neighbors:\n";
                }
            pts = annAllocPts(msh.getNbNodes(), Nodes::DIM);
            if (!pts)
                {
                std::cout << "ANN memory error while allocating points" << std::endl;
                SYSTEM_ERROR;
                }
            else if (mySets.verbose)
                {
                std::cout << "  points allocated\n";
                }
 
            for (int i = 0; i < msh.getNbNodes(); i++)
                {
                this->pts[i][0] = msh.getNode_p(i).x();
                this->pts[i][1] = msh.getNode_p(i).y();
                this->pts[i][2] = msh.getNode_p(i).z();
                }
 
            if (mySets.verbose)
                {
                std::cout << "  building kd_tree\n";
                }
 
            kdtree = new ANNkd_tree(pts, msh.getNbNodes(), Nodes::DIM);
            if (!kdtree)
                {
                std::cout << "ANN memory error while allocating kd_tree" << std::endl;
                SYSTEM_ERROR;
                }
            recenter_mem = true;
 
            if (mySets.verbose)
                {
                std::cout << "  kd_tree allocated\n";
                }
            }
        else if (mySets.verbose)
            {
            std::cout << "No recentering.\n";
            }
 
        if (mySets.restoreFileName == "")
            {
            msh.init_distrib(mySets);
            }
        else
            {
            t_prm.set_t(msh.readSol(mySets.verbose, mySets.restoreFileName));
            }
 
        /* This potentially overrides the initial time set above by msh.readSol(). */
        if (!isnan(mySets.initial_time))
            {
            t_prm.set_t(mySets.initial_time);
            }
 
        if (mySets.recenter)
            {
            direction(Nodes::IDX_Z);
            } /* DW propagation direction for recentering */
        }
 
    ~Fem()
        {
        if (recenter_mem)
            {
            annDeallocPts(pts);
            delete kdtree;
            }
        }
 
    double vmax;
 
    std::array<double,NB_ENERGY_TERMS> E;
 
    double DW_dir;
 
    double Etot0;
    
    double Etot;
 
    Mesh::mesh msh;
 
    void energy(double const t 
                ,Settings &settings );
 
    inline void evolution(void)
        {
        msh.evolution();
        Etot0 = Etot;
        }
 
    void saver(Settings &settings , timing const &t_prm ,
               std::ofstream &fout , const int nt ,
               std::vector<Eigen::Vector3d> &s  ) const;
 
    bool recenter(double thres ,
                  char recentering_direction );
 
    int time_integration(Settings &settings ,
                         LinAlgebra &linAlg ,
                         spinAcc &spinAcc_solver ,
                         scal_fmm::fmm &myFMM ,
                         timing &t_prm ,
                         int &nt );
 
private:
    bool recenter_mem;  
    ANNkd_tree *kdtree; 
    ANNpointArray pts;  
    void direction(enum Nodes::index idx_dir );
 
    void compute_all(Settings &settings , spinAcc &spinAcc_solver ,
                     scal_fmm::fmm &myFMM , const double t )
        {
        bool success(false);
        chronometer fmm_counter(2);
        myFMM.calc_demag(msh);
        if (settings.spin_acc && spinAcc_solver.compute())
            { success = true; }
        if (settings.verbose)
            {
            std::cout << "magnetostatics done in " << fmm_counter.millis() << std::endl;
            if (settings.spin_acc)
                {
                if (success)
                    { std::cout << "spin diffusion solved.\n"; }
                else
                    { std::cout << "spin diffusion solver failed.\n"; }
                }
            }
        energy(t, settings);
        evolution();
        }
    };  // end class
 
#endif