Feellgood
tetra.h
Go to the documentation of this file.
1 #ifndef tetra_h
2 #define tetra_h
3 
9 #include <set>
10 
11 #include "spinTransferTorque.h"
12 #include "facette.h"
13 #include "node.h"
14 #include "time_integration.h"
15 #include "element.h"
16 
20 namespace Tetra
21  {
22 const int N = 4;
23 const int NPI = 5;
25 const double epsilon =
26  EPSILON;
29 constexpr double A = 1. / 4.;
30 constexpr double B = 1. / 6.;
31 constexpr double C = 1. / 2.;
32 constexpr double D = -2. / 15.;
33 constexpr double E = 3. / 40.;
34 constexpr double u[NPI] = {A, B, B, B, C};
35 constexpr double v[NPI] = {A, B, B, C, B};
36 constexpr double w[NPI] = {A, B, C, B, B};
37 constexpr double pds[NPI] = {D, E, E, E, E};
45 constexpr double a[N][NPI] = {{1. - u[0] - v[0] - w[0], 1. - u[1] - v[1] - w[1],
46  1. - u[2] - v[2] - w[2], 1. - u[3] - v[3] - w[3],
47  1. - u[4] - v[4] - w[4]},
48  {u[0], u[1], u[2], u[3], u[4]},
49  {v[0], v[1], v[2], v[3], v[4]},
50  {w[0], w[1], w[2], w[3], w[4]}};
51 
53 const Eigen::Matrix<double,N,NPI> eigen_a = (Eigen::MatrixXd(N,NPI) << a[0][0], a[0][1], a[0][2], a[0][3], a[0][4],
54  a[1][0], a[1][1], a[1][2], a[1][3], a[1][4],
55  a[2][0], a[2][1], a[2][2], a[2][3], a[2][4],
56  a[3][0], a[3][1], a[3][2], a[3][3], a[3][4] ).finished();
57 
61 struct prm
62  {
63  std::string regName;
64  double alpha_LLG;
65  double A;
66  double J;
67  double K;
68  Eigen::Vector3d uk;
70  double K3;
71  Eigen::Vector3d ex;
72  Eigen::Vector3d ey;
73  Eigen::Vector3d ez;
76  double volume = 0;
79  inline void infos()
80  {
81  std::cout << "volume region name = " << regName << std::endl;
82  std::cout << "alpha_LLG = " << alpha_LLG << std::endl;
83  std::cout << "A = " << A << std::endl;
84  std::cout << "J = " << J << std::endl;
85 
86  if (K != 0)
87  {
88  std::cout << "K*uk =" << K << "*[ " << uk << "]" << std::endl;
89  }
90 
91  if (K3 != 0)
92  {
93  std::cout << "K3 = " << K3 << "; ex=[ " << ex << "], ey=[" << ey << "], ez=[" << ez
94  << "]\n";
95  }
96  };
97  };
98 
123 class Tet : public element<N,NPI>
124  {
125 public:
130  inline Tet(const std::vector<Nodes::Node> &_p_node ,
131  const int _idx ,
132  std::initializer_list<int> _i )
133  : element<N,NPI>(_p_node,_idx,_i), idx(0)
134  {
135  zeroBasing();
136  da.setZero();
137 
138  if (existNodes())
139  {
140  orientate();
141 
142  Eigen::Matrix3d J;
143  // we have to rebuild the jacobian in case of ill oriented tetrahedron
144  double detJ = Jacobian(J);
145 
146  if (fabs(detJ) < Tetra::epsilon)
147  {
148  std::cerr << "Singular jacobian in tetrahedron: |det(J)|= " << fabs(detJ) << std::endl;
149  element::infos();
150  SYSTEM_ERROR;
151  }
152  Eigen::Matrix<double,N,Nodes::DIM> dadu;
153  dadu << -1., -1., -1., 1., 0., 0., 0., 1., 0., 0., 0., 1.;
154  da = dadu * J.inverse();
155 
156  for (int j = 0; j < NPI; j++)
157  { weight[j] = detJ * Tetra::pds[j]; }
158  }
159  // do nothing lambda's (usefull for spin transfer torque)
160  extraField = [] ( Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,Tetra::NPI>> ) {};
161  extraCoeffs_BE = [](double, Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>>,
162  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>>,
163  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>>,
164  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>>,
165  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,N>>) {};
166  }
167 
169  Eigen::Matrix<double,N,Nodes::DIM> da;
170 
173  inline void interpolation(std::function<double(Nodes::Node)> getter,
174  Eigen::Ref<Eigen::Matrix<double,NPI,1>> result) const
175  {
176  Eigen::Matrix<double,N,1> scalar_nod;
177  for (int i = 0; i < N; i++) scalar_nod(i) = getter(getNode(i));
178  result = scalar_nod.transpose() * eigen_a;
179  }
180 
183  inline void interpolation(std::function<Eigen::Vector3d(Nodes::Node)> getter,
184  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> result,
185  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> Tx,
186  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> Ty,
187  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> Tz) const
188  {
189  Eigen::Matrix<double,Nodes::DIM,N> vec_nod;
190  for (int i = 0; i < N; i++) vec_nod.col(i) = getter(getNode(i));
191 
192  result = vec_nod * eigen_a;
193  Tx = vec_nod * (da.col(Nodes::IDX_X)).replicate(1,NPI);
194  Ty = vec_nod * (da.col(Nodes::IDX_Y)).replicate(1,NPI);
195  Tz = vec_nod * (da.col(Nodes::IDX_Z)).replicate(1,NPI);
196  }
197 
200  inline void interpolation_field(std::function<double(Nodes::Node)> getter,
201  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> X) const
202  {
203  Eigen::Matrix<double,N,1> scalar_nod;
204  for (int i = 0; i < N; i++) scalar_nod(i) = getter(getNode(i));
205 
206  X.setZero();
207  for (int j = 0; j < NPI; j++)
208  {
209  for (int i = 0; i < N; i++)
210  {
211  X.col(j) -= (scalar_nod[i] * da.row(i));
212  }
213  }
214  }
215 
218  inline void interpolation(std::function<double(Nodes::Node, Nodes::index)> getter, Nodes::index idx,
219  Eigen::Ref<Eigen::Matrix<double,Tetra::NPI,1>> result) const
220  {
221  Eigen::Matrix<double,N,1> scalar_nod;
222 
223  for (int i = 0; i < N; i++)
224  {
225  scalar_nod[i] = getter(getNode(i),idx);
226  }
227  result = scalar_nod.transpose() * eigen_a;
228  }
229 
239  void lumping(Eigen::Ref<Eigen::Matrix<double,NPI,1>> alpha_eff, double prefactor,
240  Eigen::Ref<Eigen::Matrix<double,3*N,3*N>> AE ) const;
241 
243  void add_drift_BE(double alpha, double s_dt, double Vdrift,
244  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> U,
245  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> V,
246  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dUd_,
247  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dVd_,
248  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,N>> BE) const;
249 
252  Eigen::Matrix<double,NPI,1> calc_aniso_uniax(Eigen::Ref<const Eigen::Vector3d> uk,
253  const double Kbis, const double s_dt,
254  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> U,
255  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> V,
256  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> H_aniso) const;
257 
260  Eigen::Matrix<double,NPI,1> calc_aniso_cub(Eigen::Ref<const Eigen::Vector3d> ex,
261  Eigen::Ref<const Eigen::Vector3d> ey,
262  Eigen::Ref<const Eigen::Vector3d> ez,
263  const double K3bis, const double s_dt,
264  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> U,
265  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> V,
266  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> H_aniso) const;
267 
271  void integrales( Tetra::prm const &param, timing const &prm_t,
272  std::function< Eigen::Matrix<double,Nodes::DIM,NPI> (void)> calc_Hext,
273  Nodes::index idx_dir, double Vdrift);
274 
276  double exchangeEnergy(Tetra::prm const &param,
277  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dudx,
278  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dudy,
279  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dudz) const;
280 
282  double anisotropyEnergy(Tetra::prm const &param, Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> u) const;
283 
285  void charges(Tetra::prm const &param,
286  std::function<Eigen::Vector3d(Nodes::Node)> getter,
287  std::vector<double> &srcDen,
288  int &nsrc) const;
289 
291  double demagEnergy(Tetra::prm const &param,
292  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dudx,
293  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dudy,
294  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dudz,
295  Eigen::Ref<Eigen::Matrix<double,NPI,1>> phi) const;
296 
298  double zeemanEnergy(Tetra::prm const &param, Eigen::Ref<Eigen::Vector3d> const Hext,
299  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> const u) const;
300 
302  double Jacobian(Eigen::Ref<Eigen::Matrix3d> J);
303 
305  double calc_vol(void) const;
306 
308  int idx;
309 
311  std::function<void( Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> H)> extraField;
312 
314  std::function<void(double Js, Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> U,
315  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dUdx,
316  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dUdy,
317  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dUdz,
318  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,N>> BE)> extraCoeffs_BE;
319 
321  void getPtGauss(Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> result) const
322  {
323  Eigen::Matrix<double,Nodes::DIM,N> vec_nod;
324  for (int i = 0; i < N; i++)
325  {
326  vec_nod.col(i) << getNode(i).p;
327  }
328  result = vec_nod * eigen_a;
329  }
330 
331 private:
333  void orientate(void)
334  {
335  if (calc_vol() < 0.0)
336  std::swap(ind[2], ind[3]);
337  }
338 
339  }; // end class Tetra
340 
343  Eigen::Matrix<double,NPI,1> calc_alpha_eff(const double dt, const double alpha,
344  Eigen::Ref<Eigen::Matrix<double,NPI,1>> uHeff);
345  } // end namespace Tetra
346 
347 #endif /* tetra_h */
Definition: tetra.h:124
void getPtGauss(Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> result) const
Definition: tetra.h:321
double calc_vol(void) const
Definition: tetra.cpp:302
void charges(Tetra::prm const &param, std::function< Eigen::Vector3d(Nodes::Node)> getter, std::vector< double > &srcDen, int &nsrc) const
Definition: tetra.cpp:245
void lumping(Eigen::Ref< Eigen::Matrix< double, NPI, 1 >> alpha_eff, double prefactor, Eigen::Ref< Eigen::Matrix< double, 3 *N, 3 *N >> AE) const
Definition: tetra.cpp:45
Eigen::Matrix< double, NPI, 1 > calc_aniso_uniax(Eigen::Ref< const Eigen::Vector3d > uk, const double Kbis, const double s_dt, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> U, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> V, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> H_aniso) const
Definition: tetra.cpp:93
std::function< void(double Js, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> U, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> dUdx, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> dUdy, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> dUdz, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, N >> BE)> extraCoeffs_BE
Definition: tetra.h:318
Eigen::Matrix< double, NPI, 1 > calc_aniso_cub(Eigen::Ref< const Eigen::Vector3d > ex, Eigen::Ref< const Eigen::Vector3d > ey, Eigen::Ref< const Eigen::Vector3d > ez, const double K3bis, const double s_dt, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> U, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> V, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> H_aniso) const
Definition: tetra.cpp:104
Eigen::Matrix< double, N, Nodes::DIM > da
Definition: tetra.h:169
Tet(const std::vector< Nodes::Node > &_p_node, const int _idx, std::initializer_list< int > _i)
Definition: tetra.h:130
double zeemanEnergy(Tetra::prm const &param, Eigen::Ref< Eigen::Vector3d > const Hext, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> const u) const
Definition: tetra.cpp:277
double exchangeEnergy(Tetra::prm const &param, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> dudx, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> dudy, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> dudz) const
Definition: tetra.cpp:212
double Jacobian(Eigen::Ref< Eigen::Matrix3d > J)
Definition: tetra.cpp:285
void add_drift_BE(double alpha, double s_dt, double Vdrift, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> U, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> V, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> dUd_, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> dVd_, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, N >> BE) const
Definition: tetra.cpp:73
int idx
Definition: tetra.h:308
std::function< void(Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> H)> extraField
Definition: tetra.h:311
void interpolation(std::function< Eigen::Vector3d(Nodes::Node)> getter, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> result, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> Tx, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> Ty, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> Tz) const
Definition: tetra.h:183
void interpolation(std::function< double(Nodes::Node, Nodes::index)> getter, Nodes::index idx, Eigen::Ref< Eigen::Matrix< double, Tetra::NPI, 1 >> result) const
Definition: tetra.h:218
void interpolation(std::function< double(Nodes::Node)> getter, Eigen::Ref< Eigen::Matrix< double, NPI, 1 >> result) const
Definition: tetra.h:173
void interpolation_field(std::function< double(Nodes::Node)> getter, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> X) const
Definition: tetra.h:200
double demagEnergy(Tetra::prm const &param, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> dudx, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> dudy, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> dudz, Eigen::Ref< Eigen::Matrix< double, NPI, 1 >> phi) const
Definition: tetra.cpp:264
void orientate(void)
Definition: tetra.h:333
double anisotropyEnergy(Tetra::prm const &param, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> u) const
Definition: tetra.cpp:223
void integrales(Tetra::prm const &param, timing const &prm_t, std::function< Eigen::Matrix< double, Nodes::DIM, NPI >(void)> calc_Hext, Nodes::index idx_dir, double Vdrift)
Definition: tetra.cpp:128
Template abstract class, mother class for tetraedrons and facettes.
Definition: element.h:25
void zeroBasing(void)
Definition: element.h:140
bool existNodes(void)
Definition: element.h:136
void infos() const
Definition: element.h:120
Eigen::Matrix< double, NPI, 1 > weight
Definition: element.h:50
std::vector< int > ind
Definition: element.h:44
const Nodes::Node & getNode(const int i) const
Definition: element.h:133
Definition: time_integration.h:6
contains namespace Facette header containing Fac class, and some constants and a less_than operator t...
index
Definition: node.h:34
const double epsilon
Definition: tetra.h:25
constexpr double A
Definition: tetra.h:29
constexpr double a[N][NPI]
Definition: tetra.h:45
constexpr double B
Definition: tetra.h:30
const int NPI
Definition: tetra.h:23
Eigen::Matrix< double, NPI, 1 > calc_alpha_eff(const double dt, const double alpha, Eigen::Ref< Eigen::Matrix< double, NPI, 1 >> uHeff)
Definition: tetra.cpp:15
const int N
Definition: tetra.h:22
constexpr double w[NPI]
Definition: tetra.h:36
constexpr double D
Definition: tetra.h:32
constexpr double C
Definition: tetra.h:31
constexpr double pds[NPI]
Definition: tetra.h:37
constexpr double v[NPI]
Definition: tetra.h:35
const Eigen::Matrix< double, N, NPI > eigen_a
Definition: tetra.h:53
constexpr double E
Definition: tetra.h:33
constexpr double u[NPI]
Definition: tetra.h:34
header to define struct Node
Definition: node.h:61
Eigen::Vector3d p
Definition: node.h:62
Definition: spinTransferTorque.h:13
Definition: tetra.h:62
Eigen::Vector3d ex
Definition: tetra.h:71
double volume
Definition: tetra.h:76
Eigen::Vector3d uk
Definition: tetra.h:68
std::string regName
Definition: tetra.h:63
double A
Definition: tetra.h:65
double alpha_LLG
Definition: tetra.h:64
Eigen::Vector3d ez
Definition: tetra.h:73
Eigen::Vector3d ey
Definition: tetra.h:72
void infos()
Definition: tetra.h:79
double J
Definition: tetra.h:66
STT p_STT
Definition: tetra.h:74
double K
Definition: tetra.h:67
double K3
Definition: tetra.h:70