-
Notifications
You must be signed in to change notification settings - Fork 4
/
tetra.h
346 lines (300 loc) · 15.7 KB
/
tetra.h
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
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
#ifndef tetra_h
#define tetra_h
/** \file tetra.h
\brief namespace Tetra
header containing Tet class, some constants, and integrales
*/
#include <set>
#include "spinTransferTorque.h"
#include "facette.h"
#include "node.h"
#include "time_integration.h"
#include "element.h"
/** \namespace Tetra
to grab altogether some constants for struct Tet
*/
namespace Tetra
{
const int N = 4; /**< number of sommits */
const int NPI = 5; /**< number of weights */
const double epsilon =
EPSILON; /**< this constant is defined from a macro in config.h.in, it is used to check the
validity of the tetrahedreon, a degeneracy test */
constexpr double A = 1. / 4.; /**< constant to build hat functions */
constexpr double B = 1. / 6.; /**< constant to build hat functions */
constexpr double C = 1. / 2.; /**< constant to build hat functions */
constexpr double D = -2. / 15.; /**< constant to build hat functions */
constexpr double E = 3. / 40.; /**< constant to build hat functions */
constexpr double u[NPI] = {A, B, B, B, C}; /**< some constants to build hat functions */
constexpr double v[NPI] = {A, B, B, C, B}; /**< some constants to build hat functions */
constexpr double w[NPI] = {A, B, C, B, B}; /**< some constants to build hat functions */
constexpr double pds[NPI] = {D, E, E, E, E}; /**< some constant weights to build hat functions */
/** constant matrix \f$ j:0..4 \f$
a[0][j] = 1.-u[j]-v[j]-w[j];
a[1][j] = u[j];
a[2][j] = v[j];
a[3][j] = w[j];
*/
constexpr double a[N][NPI] = {{1. - u[0] - v[0] - w[0], 1. - u[1] - v[1] - w[1],
1. - u[2] - v[2] - w[2], 1. - u[3] - v[3] - w[3],
1. - u[4] - v[4] - w[4]},
{u[0], u[1], u[2], u[3], u[4]},
{v[0], v[1], v[2], v[3], v[4]},
{w[0], w[1], w[2], w[3], w[4]}};
/** eigen constant matrix a */
static const Eigen::Matrix<double,N,NPI> eigen_a = [] {
Eigen::Matrix<double,N,NPI> tmp; tmp << a[0][0], a[0][1], a[0][2], a[0][3], a[0][4],
a[1][0], a[1][1], a[1][2], a[1][3], a[1][4],
a[2][0], a[2][1], a[2][2], a[2][3], a[2][4],
a[3][0], a[3][1], a[3][2], a[3][3], a[3][4];
return tmp; }();
/** \class prm
region number and material constants
*/
struct prm
{
std::string regName; /**< region name */
double alpha_LLG; /**< \f$ \alpha \f$ damping parameter */
double A; /**< exchange constant stiffness */
double J; /**< \f$ M_s = \nu_0 J \f$ */
double K; /**< uniaxial anisotropy constant */
Eigen::Vector3d uk; /**< uniaxial anisotropy axis */
double K3; /**< cubic anisotropy constant */
Eigen::Vector3d ex; /**< unit vector1 (for cubic anisotropy) */
Eigen::Vector3d ey; /**< unit vector2 (for cubic anisotropy) */
Eigen::Vector3d ez; /**< unit vector3 (for cubic anisotropy) */
STT p_STT; /**< spin transfert torque (thiaville STT) parameters */
/** print the struct parameters */
inline void infos()
{
std::cout << "volume region name = " << regName << std::endl;
std::cout << "alpha_LLG = " << alpha_LLG << std::endl;
std::cout << "A = " << A << std::endl;
std::cout << "J = " << J << std::endl;
if (K != 0)
{
std::cout << "K*uk =" << K << "*[ " << uk << "]" << std::endl;
}
if (K3 != 0)
{
std::cout << "K3 = " << K3 << "; ex=[ " << ex << "], ey=[" << ey << "], ez=[" << ez
<< "]\n";
}
};
};
/** \class Tet
Tet is a tetrahedron, containing the index references to nodes, must not be flat <br>
indices convention is<br>
```
v
.
,/
/
2(ic) 2
,/|`\ ,/|`\
,/ | `\ ,/ | `\
,/ '. `\ ,6 '. `5
,/ | `\ ,/ 8 `\
,/ | `\ ,/ | `\
0(ia)-------'.--------1(ib) --> u 0--------4--'.--------1
`\. | ,/ `\. | ,/
`\. | ,/ `\. | ,9
`\. '. ,/ `7. '. ,/
`\. |/ `\. |/
`3(id) `3
`\.
` w
```
*/
class Tet : public element<N,NPI>
{
public:
/** constructor. It initializes weight hat function with weights \f$ w_i = |J| p_i \f$
with \f$ p_i = pds[i] = (D,E,E,E,E) \f$ and dad(x|y|z) if \f$ | detJ | < \epsilon \f$ jacobian
is considered degenerated unit tests : Tet_constructor; Tet_inner_tables
*/
inline Tet(const std::vector<Nodes::Node> &_p_node /**< vector of nodes */,
const int _idx /**< [in] region index in region vector */,
std::initializer_list<int> _i /**< [in] node index */)
: element<N,NPI>(_p_node,_idx,_i), idx(0)
{
zeroBasing();
da.setZero();
if (existNodes())
{
orientate();
Eigen::Matrix3d J;
// we have to rebuild the jacobian in case of ill oriented tetrahedron
double detJ = Jacobian(J);
if (fabs(detJ) < Tetra::epsilon)
{
std::cerr << "Singular jacobian in tetrahedron: |det(J)|= " << fabs(detJ) << std::endl;
element::infos();
SYSTEM_ERROR;
}
Eigen::Matrix<double,N,Nodes::DIM> dadu;
dadu << -1., -1., -1., 1., 0., 0., 0., 1., 0., 0., 0., 1.;
da = dadu * J.inverse();
for (int j = 0; j < NPI; j++)
{ weight[j] = detJ * Tetra::pds[j]; }
}
// do nothing lambda's (usefull for spin transfer torque)
extraField = [] ( Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,Tetra::NPI>> ) {};
extraCoeffs_BE = [](double, Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>>,
Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>>,
Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>>,
Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>>,
Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,N>>) {};
}
/** local hat functions matrix, initialized by constructor: da = dadu * inverse(Jacobian) */
Eigen::Matrix<double,N,Nodes::DIM> da;
/** interpolation for scalar field : the getter function is given as a parameter in order to
* know what part of the node you want to interpolate */
inline void interpolation(std::function<double(Nodes::Node)> getter,
Eigen::Ref<Eigen::Matrix<double,NPI,1>> result) const
{
Eigen::Matrix<double,N,1> scalar_nod;
for (int i = 0; i < N; i++) scalar_nod(i) = getter(getNode(i));
result = scalar_nod.transpose() * eigen_a;
}
/** interpolation for 3D vector field and a tensor : getter function is given as a parameter to
* know what part of the node you want to interpolate */
inline 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
{
Eigen::Matrix<double,Nodes::DIM,N> vec_nod;
for (int i = 0; i < N; i++) vec_nod.col(i) = getter(getNode(i));
result = vec_nod * eigen_a;
Tx = vec_nod * (da.col(Nodes::IDX_X)).replicate(1,NPI);
Ty = vec_nod * (da.col(Nodes::IDX_Y)).replicate(1,NPI);
Tz = vec_nod * (da.col(Nodes::IDX_Z)).replicate(1,NPI);
}
/** interpolation for components of a field : the getter function is given as a parameter in
* order to know what part of the node you want to interpolate */
inline void interpolation_field(std::function<double(Nodes::Node)> getter,
Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> X) const
{
Eigen::Matrix<double,N,1> scalar_nod;
for (int i = 0; i < N; i++) scalar_nod(i) = getter(getNode(i));
X.setZero();
for (int j = 0; j < NPI; j++)
{
for (int i = 0; i < N; i++)
{
X.col(j) -= (scalar_nod[i] * da.row(i));
}
}
}
/** interpolation for the component idx of a field : the getter function is given as a parameter in order to
* know what part of the node you want to interpolate */
inline void interpolation(std::function<double(Nodes::Node, Nodes::index)> getter, Nodes::index idx,
Eigen::Ref<Eigen::Matrix<double,Tetra::NPI,1>> result) const
{
Eigen::Matrix<double,N,1> scalar_nod;
for (int i = 0; i < N; i++)
{
scalar_nod[i] = getter(getNode(i),idx);
}
result = scalar_nod.transpose() * eigen_a;
}
/** AE matrix filling with exchange contributions, diagonal contributions and magnetization contribution
AE has a block structure:
( E -Z +Y)
( +Z E -X)
( -Y +X E)
E X Y Z are N*N matrices
E is the sum of the exchange contribution and modified alpha (scheme stabilizer) on its diagonal
X Y Z are diagonal matrices with magnetization component contributions
*/
void lumping(Eigen::Ref<Eigen::Matrix<double,NPI,1>> alpha_eff, double prefactor,
Eigen::Ref<Eigen::Matrix<double,3*N,3*N>> AE ) const;
/** add drift contribution due to eventual recentering to vectors BE */
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;
/** append(+=) H_aniso for uniaxial anisotropy contribution, returns contribution to uHeff (used to
* compute the stabilizing effective damping) */
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;
/** append(+=) H_aniso for cubic anisotropy contribution, returns contribution to uHeff (used to
* compute the stabilizing effective damping) */
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;
/** computes the integral contribution of the tetrahedron to the evolution of the magnetization
*/
void integrales( Tetra::prm const ¶m, timing const &prm_t,
std::function<void( Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> Hext )> calc_Hext,
Nodes::index idx_dir, double Vdrift);
/** exchange energy of the tetrahedron */
double exchangeEnergy(Tetra::prm const ¶m,
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;
/** anisotropy energy of the tetrahedron */
double anisotropyEnergy(Tetra::prm const ¶m, Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> u) const;
/** computes volume charges, the divergence of the magnetization. Result is stored in srcDen, at nsrc position */
void charges(Tetra::prm const ¶m,
std::function<Eigen::Vector3d(Nodes::Node)> getter,
std::vector<double> &srcDen,
int &nsrc) const;
/** demagnetizing energy of the tetrahedron */
double demagEnergy(Tetra::prm const ¶m,
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;
/** zeeman energy of the tetrahedron */
double zeemanEnergy(Tetra::prm const ¶m, Eigen::Ref<Eigen::Vector3d> const Hext,
Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> const u) const;
/** \return \f$ |J| \f$ build Jacobian \f$ J \f$ */
double Jacobian(Eigen::Ref<Eigen::Matrix3d> J);
/** computes volume of the tetrahedron ; unit test Tet_calc_vol */
double calc_vol(void) const;
/** idx is the index of the tetrahedron in the vector of tetrahedron */
int idx;
/** for extra contribution to the effective field, such as spin transfert torque Hm */
std::function<void( Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> H)> extraField;
/** for extra contribution to the matrix BE, such as spin transfer torque contribs */
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;
/** returns gauss points in result = vec_nod*Tetra::a */
void getPtGauss(Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> result) const
{
Eigen::Matrix<double,Nodes::DIM,N> vec_nod;
for (int i = 0; i < N; i++)
{
vec_nod.col(i) << getNode(i).p;
}
result = vec_nod * eigen_a;
}
private:
/** orientation redefined through index swaping if needed */
void orientate(void)
{
if (calc_vol() < 0.0)
std::swap(ind[2], ind[3]);
}
}; // end class Tetra
/** to perform some second order corrections, an effective \f$ \alpha \f$ is computed here with
* a piecewise formula */
Eigen::Matrix<double,NPI,1> calc_alpha_eff(const double dt, const double alpha,
Eigen::Ref<Eigen::Matrix<double,NPI,1>> uHeff);
} // end namespace Tetra
#endif /* tetra_h */