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 static const Eigen::Matrix<double,N,NPI> eigen_a = [] {
54  Eigen::Matrix<double,N,NPI> tmp; tmp << a[0][0], a[0][1], a[0][2], a[0][3], a[0][4],
55  a[1][0], a[1][1], a[1][2], a[1][3], a[1][4],
56  a[2][0], a[2][1], a[2][2], a[2][3], a[2][4],
57  a[3][0], a[3][1], a[3][2], a[3][3], a[3][4];
58  return tmp; }();
59 
63 struct prm
64  {
65  std::string regName;
66  double alpha_LLG;
67  double A;
68  double J;
69  double K;
70  Eigen::Vector3d uk;
72  double K3;
73  Eigen::Vector3d ex;
74  Eigen::Vector3d ey;
75  Eigen::Vector3d ez;
78  double volume = 0;
81  inline void infos()
82  {
83  std::cout << "volume region name = " << regName << std::endl;
84  std::cout << "alpha_LLG = " << alpha_LLG << std::endl;
85  std::cout << "A = " << A << std::endl;
86  std::cout << "J = " << J << std::endl;
87 
88  if (K != 0)
89  {
90  std::cout << "K*uk =" << K << "*[ " << uk << "]" << std::endl;
91  }
92 
93  if (K3 != 0)
94  {
95  std::cout << "K3 = " << K3 << "; ex=[ " << ex << "], ey=[" << ey << "], ez=[" << ez
96  << "]\n";
97  }
98  };
99  };
100 
125 class Tet : public element<N,NPI>
126  {
127 public:
132  inline Tet(const std::vector<Nodes::Node> &_p_node ,
133  const int _idx ,
134  std::initializer_list<int> _i )
135  : element<N,NPI>(_p_node,_idx,_i), idx(0)
136  {
137  zeroBasing();
138  da.setZero();
139 
140  if (existNodes())
141  {
142  orientate();
143 
144  Eigen::Matrix3d J;
145  // we have to rebuild the jacobian in case of ill oriented tetrahedron
146  double detJ = Jacobian(J);
147 
148  if (fabs(detJ) < Tetra::epsilon)
149  {
150  std::cerr << "Singular jacobian in tetrahedron: |det(J)|= " << fabs(detJ) << std::endl;
151  element::infos();
152  SYSTEM_ERROR;
153  }
154  Eigen::Matrix<double,N,Nodes::DIM> dadu;
155  dadu << -1., -1., -1., 1., 0., 0., 0., 1., 0., 0., 0., 1.;
156  da = dadu * J.inverse();
157 
158  for (int j = 0; j < NPI; j++)
159  { weight[j] = detJ * Tetra::pds[j]; }
160  }
161  // do nothing lambda's (usefull for spin transfer torque)
162  extraField = [] ( Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,Tetra::NPI>> ) {};
163  extraCoeffs_BE = [](double, 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,NPI>>,
166  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>>,
167  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,N>>) {};
168  }
169 
171  Eigen::Matrix<double,N,Nodes::DIM> da;
172 
175  inline void interpolation(std::function<double(Nodes::Node)> getter,
176  Eigen::Ref<Eigen::Matrix<double,NPI,1>> result) const
177  {
178  Eigen::Matrix<double,N,1> scalar_nod;
179  for (int i = 0; i < N; i++) scalar_nod(i) = getter(getNode(i));
180  result = scalar_nod.transpose() * eigen_a;
181  }
182 
185  inline void interpolation(std::function<Eigen::Vector3d(Nodes::Node)> getter,
186  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> result,
187  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> Tx,
188  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> Ty,
189  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> Tz) const
190  {
191  Eigen::Matrix<double,Nodes::DIM,N> vec_nod;
192  for (int i = 0; i < N; i++) vec_nod.col(i) = getter(getNode(i));
193 
194  result = vec_nod * eigen_a;
195  Tx = vec_nod * (da.col(Nodes::IDX_X)).replicate(1,NPI);
196  Ty = vec_nod * (da.col(Nodes::IDX_Y)).replicate(1,NPI);
197  Tz = vec_nod * (da.col(Nodes::IDX_Z)).replicate(1,NPI);
198  }
199 
202  inline void interpolation_field(std::function<double(Nodes::Node)> getter,
203  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> X) const
204  {
205  Eigen::Matrix<double,N,1> scalar_nod;
206  for (int i = 0; i < N; i++) scalar_nod(i) = getter(getNode(i));
207 
208  X.setZero();
209  for (int j = 0; j < NPI; j++)
210  {
211  for (int i = 0; i < N; i++)
212  {
213  X.col(j) -= (scalar_nod[i] * da.row(i));
214  }
215  }
216  }
217 
220  inline void interpolation(std::function<double(Nodes::Node, Nodes::index)> getter, Nodes::index idx,
221  Eigen::Ref<Eigen::Matrix<double,Tetra::NPI,1>> result) const
222  {
223  Eigen::Matrix<double,N,1> scalar_nod;
224 
225  for (int i = 0; i < N; i++)
226  {
227  scalar_nod[i] = getter(getNode(i),idx);
228  }
229  result = scalar_nod.transpose() * eigen_a;
230  }
231 
241  void lumping(Eigen::Ref<Eigen::Matrix<double,NPI,1>> alpha_eff, double prefactor,
242  Eigen::Ref<Eigen::Matrix<double,3*N,3*N>> AE ) const;
243 
245  void add_drift_BE(double alpha, double s_dt, double Vdrift,
246  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> U,
247  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> V,
248  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dUd_,
249  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dVd_,
250  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,N>> BE) const;
251 
254  Eigen::Matrix<double,NPI,1> calc_aniso_uniax(Eigen::Ref<const Eigen::Vector3d> uk,
255  const double Kbis, const double s_dt,
256  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> U,
257  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> V,
258  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> H_aniso) const;
259 
262  Eigen::Matrix<double,NPI,1> calc_aniso_cub(Eigen::Ref<const Eigen::Vector3d> ex,
263  Eigen::Ref<const Eigen::Vector3d> ey,
264  Eigen::Ref<const Eigen::Vector3d> ez,
265  const double K3bis, const double s_dt,
266  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> U,
267  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> V,
268  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> H_aniso) const;
269 
273  void integrales( Tetra::prm const &param, timing const &prm_t,
274  std::function< Eigen::Matrix<double,Nodes::DIM,NPI> (void)> calc_Hext,
275  Nodes::index idx_dir, double Vdrift);
276 
278  double exchangeEnergy(Tetra::prm const &param,
279  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dudx,
280  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dudy,
281  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dudz) const;
282 
284  double anisotropyEnergy(Tetra::prm const &param, Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> u) const;
285 
287  void charges(Tetra::prm const &param,
288  std::function<Eigen::Vector3d(Nodes::Node)> getter,
289  std::vector<double> &srcDen,
290  int &nsrc) const;
291 
293  double demagEnergy(Tetra::prm const &param,
294  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dudx,
295  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dudy,
296  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dudz,
297  Eigen::Ref<Eigen::Matrix<double,NPI,1>> phi) const;
298 
300  double zeemanEnergy(Tetra::prm const &param, Eigen::Ref<Eigen::Vector3d> const Hext,
301  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> const u) const;
302 
304  double Jacobian(Eigen::Ref<Eigen::Matrix3d> J);
305 
307  double calc_vol(void) const;
308 
310  int idx;
311 
313  std::function<void( Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> H)> extraField;
314 
316  std::function<void(double Js, Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> U,
317  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dUdx,
318  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dUdy,
319  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dUdz,
320  Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,N>> BE)> extraCoeffs_BE;
321 
323  void getPtGauss(Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> result) const
324  {
325  Eigen::Matrix<double,Nodes::DIM,N> vec_nod;
326  for (int i = 0; i < N; i++)
327  {
328  vec_nod.col(i) << getNode(i).p;
329  }
330  result = vec_nod * eigen_a;
331  }
332 
333 private:
335  void orientate(void)
336  {
337  if (calc_vol() < 0.0)
338  std::swap(ind[2], ind[3]);
339  }
340 
341  }; // end class Tetra
342 
345  Eigen::Matrix<double,NPI,1> calc_alpha_eff(const double dt, const double alpha,
346  Eigen::Ref<Eigen::Matrix<double,NPI,1>> uHeff);
347  } // end namespace Tetra
348 
349 #endif /* tetra_h */
Definition: tetra.h:126
void getPtGauss(Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> result) const
Definition: tetra.h:323
double calc_vol(void) const
Definition: tetra.cpp:310
void charges(Tetra::prm const &param, std::function< Eigen::Vector3d(Nodes::Node)> getter, std::vector< double > &srcDen, int &nsrc) const
Definition: tetra.cpp:253
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:320
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:171
Tet(const std::vector< Nodes::Node > &_p_node, const int _idx, std::initializer_list< int > _i)
Definition: tetra.h:132
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:285
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:220
double Jacobian(Eigen::Ref< Eigen::Matrix3d > J)
Definition: tetra.cpp:293
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:310
std::function< void(Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> H)> extraField
Definition: tetra.h:313
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:185
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:220
void interpolation(std::function< double(Nodes::Node)> getter, Eigen::Ref< Eigen::Matrix< double, NPI, 1 >> result) const
Definition: tetra.h:175
void interpolation_field(std::function< double(Nodes::Node)> getter, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> X) const
Definition: tetra.h:202
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:272
void orientate(void)
Definition: tetra.h:335
double anisotropyEnergy(Tetra::prm const &param, Eigen::Ref< Eigen::Matrix< double, Nodes::DIM, NPI >> u) const
Definition: tetra.cpp:231
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
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:64
Eigen::Vector3d ex
Definition: tetra.h:73
double volume
Definition: tetra.h:78
Eigen::Vector3d uk
Definition: tetra.h:70
std::string regName
Definition: tetra.h:65
double A
Definition: tetra.h:67
double alpha_LLG
Definition: tetra.h:66
Eigen::Vector3d ez
Definition: tetra.h:75
Eigen::Vector3d ey
Definition: tetra.h:74
void infos()
Definition: tetra.h:81
double J
Definition: tetra.h:68
STT p_STT
Definition: tetra.h:76
double K
Definition: tetra.h:69
double K3
Definition: tetra.h:72