FeeLLGood/linear_algebra.cpp at master · feellgood/FeeLLGood · GitHub

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
#include "linear_algebra.h"

void LinAlgebra::base_projection()
    {
    std::mt19937 gen(rand());  // random number generator: standard Mersenne twister initialized
                               // with pseudo-random seed
    std::uniform_real_distribution<> distrib(0.0, 1.0);
    double r = distrib(gen);

    msh->setBasis(M_2_PI * r);
    }

void LinAlgebra::buildInitGuess(std::vector<double> &G) const
    {
    std::fill(G.begin(),G.end(),0);
    for (int i = 0; i < NOD; i++)
        {
        if (msh->magNode[i])
            {
            G[2*i] = msh->getProj_ep(i)/gamma0;
            G[2*i+1] = msh->getProj_eq(i)/gamma0;
            }
        }
    }

void LinAlgebra::prepareElements(Eigen::Vector3d const &Hext /**< [in] applied field */,
                                 timing const &t_prm /**< [in] */)
    {
    // it might be more efficient to build calc_Hext inside lambda of the for_each
    std::function< Eigen::Matrix<double,Nodes::DIM,Tetra::NPI>(void)> calc_Hext =
                  [&Hext](void)
                      {
                      Eigen::Matrix<double,Nodes::DIM,Tetra::NPI> H;
                      H.colwise() = Hext;  // set all columns of H to Hext
                      return H;
                      };

    std::for_each(EXEC_POL, msh->tet.begin(), msh->tet.end(),
                  [this, &calc_Hext, &t_prm](Tetra::Tet &tet)
                  {
                  if (msh->isMagnetic(tet))
                      { tet.integrales(paramTet[tet.idxPrm], t_prm, calc_Hext, idx_dir, DW_vz); }
                  });

    std::for_each(EXEC_POL, msh->fac.begin(), msh->fac.end(),
                  [this](Facette::Fac &fac)
                  {
                  // the contribution to Lp computed in integrales is due to surface anisotropy
                  if ((msh->isMagnetic(fac))&&(paramFac[fac.idxPrm].Ks != 0))
                      { fac.integrales(paramFac[fac.idxPrm]);  }
                  });
    }

void LinAlgebra::prepareElements(double const A_Hext /**< [in] amplitude applied field (might be time dependant)*/,
                                 timing const &t_prm /**< [in] */)
    {
    std::for_each(EXEC_POL, msh->tet.begin(), msh->tet.end(),
                  [this, &A_Hext, &t_prm](Tetra::Tet &tet)
                  {
                  if (msh->isMagnetic(tet))
                      {
                      Eigen::Matrix<double,Nodes::DIM,Tetra::NPI> sp_H = extSpaceField[tet.idx];
                      std::function< Eigen::Matrix<double,Nodes::DIM,Tetra::NPI> (void)> calc_Hext =
                      [&sp_H,&A_Hext](void)
                          {
                          Eigen::Matrix<double,Nodes::DIM,Tetra::NPI> H= A_Hext*sp_H;
                          return H;
                          };
                      tet.integrales(paramTet[tet.idxPrm], t_prm, calc_Hext, idx_dir, DW_vz);
                      }
                  });

    std::for_each(EXEC_POL, msh->fac.begin(), msh->fac.end(),
                  [this](Facette::Fac &fac)
                  {
                  // the contribution to Lp computed in integrales is due to surface anisotropy
                  if ((msh->isMagnetic(fac))&&(paramFac[fac.idxPrm].Ks != 0))
                    { fac.integrales(paramFac[fac.idxPrm]);  }
                  });
    }