v0.14.0
Loading...
Searching...
No Matches
ElasticOperators.hpp
Go to the documentation of this file.
1/* This file is part of MoFEM.
2 * MoFEM is free software: you can redistribute it and/or modify it under
3 * the terms of the GNU Lesser General Public License as published by the
4 * Free Software Foundation, either version 3 of the License, or (at your
5 * option) any later version.
6 *
7 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
8 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
9 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
10 * License for more details.
11 *
12 * You should have received a copy of the GNU Lesser General Public
13 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
14
15namespace OpElasticTools {
16
17// constexpr auto VecSetValues = MoFEM::VecSetValues<EssentialBcStorage>;
18// constexpr auto MatSetValues = MoFEM::MatSetValues<EssentialBcStorage>;
19// constexpr auto MatSetValues = MoFEM::MatSetValues<MpEntityStorage>;
20
21using OpMassPrecHTensorBound = FormsIntegrators<BoundaryEleOp>::Assembly<
22 USER_ASSEMBLE>::BiLinearForm<GAUSS>::OpMass<3, 9>;
23
24using OpMassPrecHTensorVol = FormsIntegrators<DomainEleOp>::Assembly<
25 USER_ASSEMBLE>::BiLinearForm<GAUSS>::OpMass<3, 9>;
26
27using OpMassPrecScalar = FormsIntegrators<DomainEleOp>::Assembly<
28 USER_ASSEMBLE>::BiLinearForm<GAUSS>::OpMass<1, 1>;
29
30//! [Common data]
31struct CommonData {
32 // Ddg<double, 3, 3> tD;
33 boost::shared_ptr<MatrixDouble> mtD;
34 boost::shared_ptr<MatrixDouble> mtD_Axiator;
35 boost::shared_ptr<MatrixDouble> mtD_Deviator;
36
37 boost::shared_ptr<MatrixDouble> mGradPtr;
38 boost::shared_ptr<MatrixDouble> mStrainPtr;
39 boost::shared_ptr<MatrixDouble> mStressPtr;
40 boost::shared_ptr<MatrixDouble> mPiolaStressPtr;
41
42 boost::shared_ptr<MatrixDouble> mDeformationGradient;
43 boost::shared_ptr<MatrixDouble> matC;
44 boost::shared_ptr<VectorDouble> detF;
45
46 boost::shared_ptr<MatrixDouble> matEigenVal;
47 boost::shared_ptr<MatrixDouble> matEigenVector;
48
49 boost::shared_ptr<MatrixDouble> materialTangent;
50 boost::shared_ptr<MatrixDouble> dE_dF_mat;
51
52 array<int, 2> stateArrayPlast;
53 array<int, 2> stateArrayCont;
54 array<double, 2> stateArrayArea;
55 vector<double> bodyReactions;
56
57 SmartPetscObj<Vec> reactionsVec;
59
61 double mMeasure;
62 double time;
63
64 VectorInt rowIndices;
65
66 boost::shared_ptr<MatrixDouble> mDispPtr;
68
70
72 : stateArrayCont({0, 0}), stateArrayPlast({0, 0}),
73 stateArrayArea({0., 0.}), isNeoHook(false) {
74 }
75
77
79
80 SmartPetscObj<Mat> SEpEp;
81 SmartPetscObj<AO> aoSEpEp;
82 SmartPetscObj<Mat> STauTau;
83 SmartPetscObj<AO> aoSTauTau;
84
85 boost::shared_ptr<BlockMatContainer> blockMatContainerPtr;
86 MatrixDouble massMatTensor;
87 MatrixDouble massMatScalar;
88};
89
90template <typename T> inline auto to_non_symm(T &symm) {
91 Tensor2<double, 3, 3> non_symm;
92 Number<0> N0;
93 Number<1> N1;
94 Number<2> N2;
95 non_symm(N0, N0) = symm(N0, N0);
96 non_symm(N1, N1) = symm(N1, N1);
97 non_symm(N2, N2) = symm(N2, N2);
98 non_symm(N0, N1) = non_symm(N1, N0) = symm(N0, N1);
99 non_symm(N0, N2) = non_symm(N2, N0) = symm(N0, N2);
100 non_symm(N2, N1) = non_symm(N1, N2) = symm(N2, N1);
101 return non_symm;
102};
103
104template <typename T> inline auto to_symm(T &non_symm) {
105 Tensor2_symmetric<double, 3> symm;
106 Number<0> N0;
107 Number<1> N1;
108 Number<2> N2;
109 symm(N0, N0) = non_symm(N0, N0);
110 symm(N1, N1) = non_symm(N1, N1);
111 symm(N2, N2) = non_symm(N2, N2);
112 symm(N0, N1) = non_symm(N0, N1);
113 symm(N0, N2) = non_symm(N0, N2);
114 symm(N1, N2) = non_symm(N1, N2);
115 return symm;
116};
117
118inline bool is_eq(const double &a, const double &b) {
119 constexpr double eps = 1e-12;
120 return abs(a - b) < eps;
121};
122inline bool is_diff(const double &a, const double &b) {
123 constexpr double eps = 1e-12;
124 return abs(a - b) > eps;
125};
126
127template <typename T>
128inline Ddg<double, 3, 3> tensor4_to_ddg(Tensor4<T, 3, 3, 3, 3> &t) {
129 Ddg<double, 3, 3> loc_tens;
130
131 Number<0> N0;
132 Number<1> N1;
133 Number<2> N2;
134
135 loc_tens(N0, N0, N0, N0) = t(N0, N0, N0, N0); // 0
136 loc_tens(N0, N0, N0, N1) = t(N0, N0, N0, N1); // 1
137 loc_tens(N0, N0, N0, N2) = t(N0, N0, N0, N2); // 2
138 loc_tens(N0, N0, N1, N1) = t(N0, N0, N1, N1); // 3
139 loc_tens(N0, N0, N1, N2) = t(N0, N0, N1, N2); // 4
140 loc_tens(N0, N0, N2, N2) = t(N0, N0, N2, N2); // 5
141 loc_tens(N0, N1, N0, N0) = t(N0, N1, N0, N0); // 6
142 loc_tens(N0, N1, N0, N1) = t(N0, N1, N0, N1); // 7
143 loc_tens(N0, N1, N0, N2) = t(N0, N1, N0, N2); // 8
144 loc_tens(N0, N1, N1, N1) = t(N0, N1, N1, N1); // 9
145 loc_tens(N0, N1, N1, N2) = t(N0, N1, N1, N2); // 10
146 loc_tens(N0, N1, N2, N2) = t(N0, N1, N2, N2); // 11
147 loc_tens(N0, N2, N0, N0) = t(N0, N2, N0, N0); // 12
148 loc_tens(N0, N2, N0, N1) = t(N0, N2, N0, N1); // 13
149 loc_tens(N0, N2, N0, N2) = t(N0, N2, N0, N2); // 14
150 loc_tens(N0, N2, N1, N1) = t(N0, N2, N1, N1); // 15
151 loc_tens(N0, N2, N1, N2) = t(N0, N2, N1, N2); // 16
152 loc_tens(N0, N2, N2, N2) = t(N0, N2, N2, N2); // 17
153 loc_tens(N1, N1, N0, N0) = t(N1, N1, N0, N0); // 18
154 loc_tens(N1, N1, N0, N1) = t(N1, N1, N0, N1); // 19
155 loc_tens(N1, N1, N0, N2) = t(N1, N1, N0, N2); // 20
156 loc_tens(N1, N1, N1, N1) = t(N1, N1, N1, N1); // 21
157 loc_tens(N1, N1, N1, N2) = t(N1, N1, N1, N2); // 22
158 loc_tens(N1, N1, N2, N2) = t(N1, N1, N2, N2); // 23
159 loc_tens(N1, N2, N0, N0) = t(N1, N2, N0, N0); // 24
160 loc_tens(N1, N2, N0, N1) = t(N1, N2, N0, N1); // 25
161 loc_tens(N1, N2, N0, N2) = t(N1, N2, N0, N2); // 26
162 loc_tens(N1, N2, N1, N1) = t(N1, N2, N1, N1); // 27
163 loc_tens(N1, N2, N1, N2) = t(N1, N2, N1, N2); // 28
164 loc_tens(N1, N2, N2, N2) = t(N1, N2, N2, N2); // 29
165 loc_tens(N2, N2, N0, N0) = t(N2, N2, N0, N0); // 30
166 loc_tens(N2, N2, N0, N1) = t(N2, N2, N0, N1); // 31
167 loc_tens(N2, N2, N0, N2) = t(N2, N2, N0, N2); // 32
168 loc_tens(N2, N2, N1, N1) = t(N2, N2, N1, N1); // 33
169 loc_tens(N2, N2, N1, N2) = t(N2, N2, N1, N2); // 34
170 loc_tens(N2, N2, N2, N2) = t(N2, N2, N2, N2); // 35
171
172 return loc_tens;
173};
174
175template <typename T> inline auto get_tensor4_from_ddg(Ddg<T, 3, 3> &ddg) {
176 Tensor4<double, 3, 3, 3, 3> tens4;
177 for (int ii = 0; ii != 3; ++ii)
178 for (int jj = 0; jj != 3; ++jj)
179 for (int kk = 0; kk != 3; ++kk)
180 for (int ll = 0; ll != 3; ++ll)
181 tens4(ii, jj, kk, ll) = ddg(ii, jj, kk, ll);
182 return tens4;
183};
184
185template <typename T> inline bool is_ddg(Tensor4<T, 3, 3, 3, 3> tensor) {
186 auto test_ddg = tensor4_to_ddg(tensor);
187 for (int ii = 0; ii != 3; ++ii)
188 for (int jj = 0; jj != 3; ++jj)
189 for (int kk = 0; kk != 3; ++kk)
190 for (int ll : {0, 1, 2}) {
191 const double diff =
192 abs(test_ddg(ii, jj, kk, ll) - tensor(ii, jj, kk, ll));
193 if (diff > 1e-6)
194 return false;
195 }
196
197 return true;
198};
199
200template <typename T> inline auto exp_map(const Tensor2_symmetric<T, 3> &X) {
201 Tensor2<double, 3, 3> exp_tensor(1., 0., 0., 0., 1., 0., 0., 0., 1.);
202 Tensor2<double, 3, 3> X_nt;
203
204 constexpr int max_it = 1000;
205 constexpr double eps = 1e-12;
206
207 int fac_n = 1;
208
209 auto Xcopy = to_non_symm(X);
210 auto X_n = to_non_symm(X);
211 exp_tensor(i, j) += Xcopy(i, j);
212
213 for (int n = 2; n != max_it; ++n) {
214 fac_n *= n;
215 X_nt(i, j) = X_n(i, j);
216 X_n(i, k) = X_nt(i, j) * Xcopy(j, k);
217 const double inv_nf = 1. / fac_n;
218 exp_tensor(i, j) += inv_nf * X_n(i, j);
219 const double e = inv_nf * std::abs(sqrt(X_n(i, j) * X_n(i, j)));
220 if (e < eps)
221 return to_symm(exp_tensor);
222 }
223
224 return to_symm(exp_tensor);
225};
226
227template <typename T>
228inline auto get_rotation_R(Tensor1<T, 3> t1_omega, double tt) {
229
230 auto get_rotation = [&](Tensor1<double, 3> &t_omega) {
232
233 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
234 t_R(i, j) = t_kd(i, j);
235
236 const double angle = sqrt(t_omega(i) * t_omega(i));
237 if (std::abs(angle) < 1e-18)
238 return t_R;
239
241 t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
242 const double a = sin(angle) / angle;
243 const double ss_2 = sin(angle / 2.);
244 const double b = 2. * ss_2 * ss_2 / (angle * angle);
245 t_R(i, j) += a * t_Omega(i, j);
246 t_R(i, j) += b * t_Omega(i, k) * t_Omega(k, j);
247
248 return t_R;
249 };
250
251 Tensor1<double, 3> my_c_omega;
252 my_c_omega(i) = t1_omega(i) * tt;
253 auto t_R = get_rotation(my_c_omega);
254
255 return get_rotation(my_c_omega);
256};
257
258template <typename T>
259inline auto get_log_map(const Tensor2_symmetric<T, 3> &X) {
260
261 constexpr int max_it = 1000;
262 constexpr double eps = 1e-12;
263
264 auto log_tensor = to_non_symm(X);
265 log_tensor(i, j) -= kronecker_delta(i, j);
266
267 auto Xcopy = to_non_symm(log_tensor);
268 auto X_n = to_non_symm(log_tensor);
269 Tensor2<double, 3, 3> X_nt_tmp;
270
271 for (int n = 2; n != max_it; ++n) {
272
273 X_nt_tmp(i, j) = X_n(i, j);
274 X_n(i, k) = X_nt_tmp(i, j) * Xcopy(j, k);
275 const double inv_n = 1. / n;
276 log_tensor(i, j) += inv_n * pow(-1, n + 1) * X_n(i, j);
277 const double e = inv_n * std::abs(sqrt(X_n(i, j) * X_n(i, j)));
278 if (e < eps)
279 break;
280 }
281
282 return to_symm(log_tensor);
283};
284
285template <typename T> inline auto get_diff_log_map(Tensor2_symmetric<T, 3> &X) {
286 Tensor2<double, 3, 3> X_nt_tmp;
287 Tensor4<double, 3, 3, 3, 3> diff_log_map;
288 Tensor4<double, 3, 3, 3, 3> diff_log_map_tmp;
289
290 constexpr int max_it = 1000;
291 constexpr double eps = 1e-12;
292
293 auto get_Xn = [&](auto &log_tensor) {
294 vector<Tensor2<double, 3, 3>> Xn;
295 Xn.emplace_back();
296 auto &cu_Xn = Xn.back();
297 cu_Xn(i, j) = kronecker_delta(i, j);
298 for (int n = 2; n != max_it; ++n) {
299 const double cof = 1. / n;
300 auto X_nt_tmp = Xn.back();
301 Xn.emplace_back();
302 auto &c_Xn = Xn.back();
303 c_Xn(i, k) = X_nt_tmp(i, j) * log_tensor(j, k);
304 const double e = cof * std::abs(sqrt(c_Xn(i, j) * c_Xn(i, j)));
305 if (e < eps)
306 return Xn;
307 }
308 return Xn;
309 };
310
311 auto log_tensor = to_non_symm(X);
312 log_tensor(i, j) -= kronecker_delta(i, j);
313 diff_log_map(i, j, k, l) = 0.;
314
315 auto Xn_vec = get_Xn(log_tensor);
316
317 for (int n = 1; n != Xn_vec.size() + 1; ++n) {
318 const double cof = pow(-1, n + 1) / n;
319 for (int m = 1; m != n + 1; ++m) {
320 auto &Xn_1 = Xn_vec[m - 1];
321 auto &Xn_2 = Xn_vec[n - m];
322 diff_log_map(i, j, k, l) += cof * Xn_1(i, k) * Xn_2(l, j);
323 }
324 }
325
326 auto isddg = tensor4_to_ddg((*cache).Is);
327 diff_log_map_tmp(i, j, k, l) = isddg(i, j, m, n) * diff_log_map(m, n, k, l);
328 return tensor4_to_ddg(diff_log_map_tmp);
329};
330
331template <typename T>
332inline MoFEMErrorCode
333get_eigen_val_and_proj_lapack(const Tensor2_symmetric<T, 3> &X,
335 Tensor2<T, 3, 3> &eig_vec) {
336
339 auto eigen_vec_lap = to_non_symm(X);
340 constexpr int n = 3;
341 constexpr int lda = 3;
342 constexpr int lwork = 15;
343 std::array<double, 15> work;
344 if (lapack_dsyev('V', 'U', n, &(eigen_vec_lap(0, 0)), lda, &eig(0),
345 work.data(), lwork) > 0)
346 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
347 "The algorithm failed to compute eigenvalues.");
348
349 lambda(i) = eig(i);
350 eig_vec(i, j) = eigen_vec_lap(i, j);
351
352 // perturb
353 // const double h = 1e-12;
354 // if (is_eq(lambda[0], lambda[1]))
355 // lambda[0] *= (1 + h);
356 // if (is_eq(lambda[1], lambda[2]))
357 // lambda[2] *= (1 - h);
358
360}
361
362template <typename T1, typename T2, typename T3>
363inline auto get_ln_X_lapack(const Tensor2_symmetric<T1, 3> &X,
365 Tensor2<T3, 3, 3> &eig_vec) {
366 Tensor2_symmetric<double, 3> lnX;
367 Tensor2_symmetric<double, 3> diag(0., 0., 0., 0., 0., 0.);
368 for (int ii = 0; ii != 3; ++ii)
369 diag(ii, ii) = log(lambda(ii));
370
371 lnX(i, j) = eig_vec(l, i) ^ (diag(l, k) * eig_vec(k, j));
372
373 return lnX;
374};
375
376template <typename T> inline auto get_ln_X(const Tensor2_symmetric<T, 3> &X) {
377 return get_log_map(X);
378};
379
380template <typename T, typename TP, typename T3>
381inline auto
382calculate_nominal_moduli_TL(Tensor2_symmetric<T, 3> &X, Tensor2<TP, 3, 3> &F,
383 Tensor2_symmetric<TP, 3> &stress,
384 Tensor1<T, 3> &lambda, Tensor2<T, 3, 3> &eig_vec,
385 Tensor4<T3, 3, 3, 3, 3> &TLs3,
386 bool compute_tangent) {
387
388 // Tensor4<double, 3, 3, 3, 3> TLs3;
389 Tensor4<double, 3, 3, 3, 3> P3;
390 Tensor2<double, 3, 3> invF;
391
392 double det_f;
393 CHKERR determinantTensor3by3(F, det_f);
394 CHKERR invertTensor3by3(F, det_f, invF);
395
396 // from https://www.sciencedirect.com/science/article/pii/S0045782502004383
397
398 // VectorDouble lambda(3);
399 VectorDouble e(3);
400 VectorDouble f(3);
401 VectorDouble d(3);
402 MatrixDouble Vi(3, 3);
403 MatrixDouble Ksi(3, 3);
404 MatrixDouble Zeta(3, 3);
405 Vi.clear();
406 Ksi.clear();
407 Zeta.clear();
408 double eta = 0;
409
410 Tensor1<double, 3> _N[3];
411 Tensor1<double, 3> _n[3];
412
413 for (int dd : {0, 1, 2}) {
414 for (int cc : {0, 1, 2})
415 _N[dd](cc) = eig_vec(dd, cc);
416 }
417
418 for (int dd : {0, 1, 2})
419 _n[dd](i) = F(i, j) * _N[dd](j);
420
421 for (int ii = 0; ii != 3; ++ii) {
422 e(ii) = 0.5 * log(lambda(ii));
423 d(ii) = 1. / lambda(ii);
424 f(ii) = -2. / (lambda(ii) * lambda(ii));
425 }
426
427 const double &lam0 = lambda(0);
428 const double &lam1 = lambda(1);
429 const double &lam2 = lambda(2);
430
431 auto compute_ksi_eta_2eq = [&](int i_, int j_, int k_) {
432 Vi(i_, j_) = Vi(j_, i_) = 0.5 * d(i_);
433 Ksi(i_, j_) = Ksi(j_, i_) = 0.125 * f(i_);
434
435 for (int ii = 0; ii != 3; ++ii)
436 for (int jj = 0; jj != 3; ++jj)
437 if (ii != jj) {
438 if (jj == k_ || ii == k_) {
439 Vi(ii, jj) = (e(ii) - e(jj)) / (lambda(ii) - lambda(jj));
440 Ksi(ii, jj) =
441 (Vi(ii, jj) - 0.5 * d(jj)) / (lambda(ii) - lambda(jj));
442 }
443 }
444 eta = Ksi(k_, i_);
445 };
446
447 if (is_diff(lam0, lam1) && is_diff(lam0, lam2) && is_diff(lam1, lam2)) {
448 // all different
449 for (int ii = 0; ii != 3; ++ii)
450 for (int jj = 0; jj != 3; ++jj) {
451 if (ii != jj) {
452
453 Vi(ii, jj) = (e(ii) - e(jj)) / (lambda(ii) - lambda(jj));
454 Ksi(ii, jj) = (Vi(ii, jj) - 0.5 * d(jj)) / (lambda(ii) - lambda(jj));
455 for (int kk = 0; kk != 3; ++kk) {
456
457 if (kk != ii && kk != jj)
458 eta += e(ii) / (2. * (lambda(ii) - lambda(jj)) *
459 (lambda(ii) - lambda(kk)));
460 }
461 }
462 }
463
464 } else if (is_eq(lam0, lam1) && is_eq(lam0, lam2) && is_eq(lam1, lam2)) {
465 // all the same
466 for (int ii = 0; ii != 3; ++ii)
467 for (int jj = 0; jj != 3; ++jj)
468 if (ii != jj) {
469
470 Vi(ii, jj) = 0.5 * d(ii);
471 Ksi(ii, jj) = 0.125 * f(ii);
472 }
473 eta = 0.125 * f(0);
474 } else if (is_eq(lam0, lam1)) {
475
476 compute_ksi_eta_2eq(0, 1, 2);
477
478 } else if (is_eq(lam0, lam2)) {
479
480 compute_ksi_eta_2eq(0, 2, 1);
481
482 } else if (is_eq(lam1, lam2)) {
483
484 compute_ksi_eta_2eq(1, 2, 0);
485 }
486
487 P3(i, j, k, l) = 0.;
488
489 Tensor2<double, 3, 3> Mii;
490 Tensor2<double, 3, 3> Mij;
491 Tensor2<double, 3, 3> Mik;
492 Tensor2<double, 3, 3> Mjk;
493 Tensor2<double, 3, 3> Mjj;
494
495 for (int ii = 0; ii != 3; ++ii) {
496 Mii(i, j) = _n[ii](i) * _N[ii](j) + _n[ii](i) * _N[ii](j);
497 P3(i, j, k, l) += 0.5 * d(ii) * _N[ii](i) * _N[ii](j) * Mii(k, l);
498 for (int jj = 0; jj != 3; ++jj)
499 if (ii != jj) {
500 Mij(i, j) = _n[ii](i) * _N[jj](j) + _n[jj](i) * _N[ii](j);
501 P3(i, j, k, l) += Vi(ii, jj) * _N[ii](i) * _N[jj](j) * Mij(k, l);
502 }
503 }
504
505 if (!compute_tangent)
506 return P3;
507
508 TLs3(i, j, k, l) = 0.;
509 for (int ii = 0; ii != 3; ++ii)
510 for (int jj = 0; jj != 3; ++jj)
511 Zeta(ii, jj) = (stress(i, j) * (_N[ii](i) * _N[jj](j)));
512
513 for (int ii = 0; ii != 3; ++ii) {
514 Mii(i, j) = _n[ii](i) * _N[ii](j) + _n[ii](i) * _N[ii](j);
515 TLs3(i, j, k, l) += 0.25 * f(ii) * Zeta(ii, ii) * Mii(i, j) * Mii(k, l);
516 }
517
518 for (int ii = 0; ii != 3; ++ii)
519 for (int jj = 0; jj != 3; ++jj)
520 if (jj != ii)
521 for (int kk = 0; kk != 3; ++kk)
522 if (kk != ii && kk != jj) {
523 Mik(i, j) = _n[ii](i) * _N[kk](j) + _n[kk](i) * _N[ii](j);
524 Mjk(i, j) = _n[jj](i) * _N[kk](j) + _n[kk](i) * _N[jj](j);
525 TLs3(i, j, k, l) += 2. * eta * Zeta(ii, jj) * Mik(i, j) * Mjk(k, l);
526 }
527
528 for (int ii = 0; ii != 3; ++ii)
529 for (int jj = 0; jj != 3; ++jj)
530 if (jj != ii) {
531 Mij(i, j) = _n[ii](i) * _N[jj](j) + _n[jj](i) * _N[ii](j);
532 Mjj(i, j) = _n[jj](i) * _N[jj](j) + _n[jj](i) * _N[jj](j);
533 const double k2 = 2. * Ksi(ii, jj);
534 const double zp = k2 * Zeta(ii, jj);
535 const double zp2 = k2 * Zeta(jj, jj);
536 TLs3(i, j, k, l) += zp * Mij(i, j) * Mjj(k, l) +
537 zp * Mjj(i, j) * Mij(k, l) +
538 zp2 * Mij(i, j) * Mij(k, l);
539 }
540
541 Tensor2<double, 3, 3> Piola;
542 Tensor2<double, 3, 3> invFPiola;
543 Tensor4<double, 3, 3, 3, 3> Ttmp3;
544 auto nstress = to_non_symm(stress);
545 Piola(k, l) = nstress(m, n) * P3(m, n, k, l);
546 invFPiola(i, j) = invF(i, k) * Piola(k, j);
547 Ttmp3(i, j, k, l) = 0.5 * invFPiola(i, k) * kronecker_delta(j, l) +
548 0.5 * invFPiola(i, l) * kronecker_delta(j, k);
549 TLs3(i, j, k, l) += Ttmp3(i, j, k, l);
550
551 return P3;
552};
553
554template <typename T, typename TP, typename T3>
555inline auto
556test_projection_lukasz(Tensor2_symmetric<T, 3> &X,
557 Tensor2_symmetric<TP, 3> &stress, Tensor1<T, 3> &lambda,
558 Tensor2<T, 3, 3> &eig_vec, Ddg<T3, 3, 3> &TLs3) {
559
560 Tensor1<double, 3> eig{lambda(0), lambda(1), lambda(2)};
561 Tensor2<double, 3, 3> eigen_vec;
562 eigen_vec(i, j) = eig_vec(i, j);
563
564 auto get_uniq_nb = [&](auto eig) {
565 return distance(&eig(0), unique(&eig(0), &eig(0) + 3, is_eq));
566 };
567
568 auto sort_eigen_vals2 = [&](auto &eig, auto &eigen_vec) {
569 if (is_eq(eig(0), eig(1))) {
570 Tensor2<double, 3, 3> eigen_vec_c{
571 eigen_vec(0, 0), eigen_vec(0, 1), eigen_vec(0, 2),
572 eigen_vec(2, 0), eigen_vec(2, 1), eigen_vec(2, 2),
573 eigen_vec(1, 0), eigen_vec(1, 1), eigen_vec(1, 2)};
574 Tensor1<double, 3> eig_c{eig(0), eig(2), eig(1)};
575 eig(i) = eig_c(i);
576 eigen_vec(i, j) = eigen_vec_c(i, j);
577 } else {
578 Tensor2<double, 3, 3> eigen_vec_c{
579 eigen_vec(1, 0), eigen_vec(1, 1), eigen_vec(1, 2),
580 eigen_vec(0, 0), eigen_vec(0, 1), eigen_vec(0, 2),
581 eigen_vec(2, 0), eigen_vec(2, 1), eigen_vec(2, 2)};
582 Tensor1<double, 3> eig_c{eig(1), eig(0), eig(2)};
583 eig(i) = eig_c(i);
584 eigen_vec(i, j) = eigen_vec_c(i, j);
585 }
586 };
587
588 Tensor2_symmetric<double, 3> Ts;
589 Ts(i, j) = stress(i, j);
590 auto f = [](double v) { return 0.5 * log(v); };
591 auto d_f = [](double v) { return 0.5 / v; };
592 auto dd_f = [](double v) { return -0.5 / (v * v); };
593
594 auto nb_uniq = get_uniq_nb(eig);
595 if (nb_uniq == 2)
596 sort_eigen_vals2(eig, eigen_vec);
597 auto t_d =
598
599 EigenMatrix::getDiffMat(eig, eigen_vec, f, d_f, nb_uniq);
600
601 auto t_dd =
602
603 EigenMatrix::getDiffDiffMat(eig, eigen_vec, f, d_f, dd_f, Ts, nb_uniq);
604
605 TLs3(i, j, k, l) = t_dd(i, j, k, l);
606 return t_d;
607};
608
609template <bool TANGENT, typename T, typename TP, typename T3>
610inline auto calculate_lagrangian_moduli_TL(Tensor2_symmetric<TP, 3> &stress,
612 Tensor2<T, 3, 3> &eig_vec,
613 Ddg<T3, 3, 3> &TLs3) {
614
615 Ddg<double, 3, 3> P3;
616 // from https://www.sciencedirect.com/science/article/pii/S0045782502004383
617
618 VectorDouble e(3);
619 VectorDouble f(3);
620 VectorDouble d(3);
621 MatrixDouble Vi(3, 3);
622 MatrixDouble Ksi(3, 3);
623 MatrixDouble Zeta(3, 3);
624 Vi.clear();
625 Ksi.clear();
626 Zeta.clear();
627 double eta = 0;
628
629 Tensor1<double, 3> _N[3];
630 for (int dd = 0; dd != 3; ++dd) {
631 for (int cc = 0; cc != 3; ++cc)
632 _N[dd](cc) = eig_vec(dd, cc);
633 }
634
635 for (int ii = 0; ii != 3; ++ii) {
636 e(ii) = 0.5 * log(lambda(ii));
637 d(ii) = 1. / lambda(ii);
638 f(ii) = -2. / (lambda(ii) * lambda(ii));
639 }
640
641 const double &lam0 = lambda(0);
642 const double &lam1 = lambda(1);
643 const double &lam2 = lambda(2);
644
645 auto compute_ksi_eta_2eq = [&](int i_, int j_, int k_) {
646 Vi(i_, j_) = Vi(j_, i_) = 0.5 * d(i_);
647 Ksi(i_, j_) = Ksi(j_, i_) = 0.125 * f(i_);
648
649 for (int ii = 0; ii != 3; ++ii)
650 for (int jj = 0; jj != 3; ++jj)
651 if (ii != jj) {
652 if (jj == k_ || ii == k_) {
653 Vi(ii, jj) = (e(ii) - e(jj)) / (lambda(ii) - lambda(jj));
654 Ksi(ii, jj) =
655 (Vi(ii, jj) - 0.5 * d(jj)) / (lambda(ii) - lambda(jj));
656 }
657 }
658 eta = Ksi(k_, i_);
659 };
660
661 if (is_diff(lam0, lam1) && is_diff(lam0, lam2) && is_diff(lam1, lam2)) {
662 // all different
663 for (int ii = 0; ii != 3; ++ii)
664 for (int jj = 0; jj != 3; ++jj) {
665 if (ii != jj) {
666
667 Vi(ii, jj) = (e(ii) - e(jj)) / (lambda(ii) - lambda(jj));
668 Ksi(ii, jj) = (Vi(ii, jj) - 0.5 * d(jj)) / (lambda(ii) - lambda(jj));
669 for (int kk = 0; kk != 3; ++kk) {
670
671 if (kk != ii && kk != jj)
672 eta += e(ii) / (2. * (lambda(ii) - lambda(jj)) *
673 (lambda(ii) - lambda(kk)));
674 }
675 }
676 }
677
678 } else if (is_eq(lam0, lam1) && is_eq(lam0, lam2) && is_eq(lam1, lam2)) {
679 // all the same
680 for (int ii = 0; ii != 3; ++ii)
681 for (int jj = 0; jj != 3; ++jj)
682 if (ii != jj) {
683
684 Vi(ii, jj) = 0.5 * d(ii);
685 Ksi(ii, jj) = 0.125 * f(ii);
686 }
687 eta = 0.125 * f(0);
688 } else if (is_eq(lam0, lam1)) {
689
690 compute_ksi_eta_2eq(0, 1, 2);
691
692 } else if (is_eq(lam0, lam2)) {
693
694 compute_ksi_eta_2eq(0, 2, 1);
695
696 } else if (is_eq(lam1, lam2)) {
697
698 compute_ksi_eta_2eq(1, 2, 0);
699 }
700
701 P3(i, j, k, l) = 0.;
702 Tensor2_symmetric<double, 3> Mii;
703 Tensor2_symmetric<double, 3> Mij;
704 Tensor2_symmetric<double, 3> Mik;
705 Tensor2_symmetric<double, 3> Mjk;
706 Tensor2_symmetric<double, 3> Mjj;
707
708 for (int ii = 0; ii != 3; ++ii) {
709 Mii(i, j) = (_N[ii](i) ^ _N[ii](j)) + (_N[ii](i) ^ _N[ii](j));
710 P3(i, j, k, l) += ((0.5 * d(ii) * _N[ii](i)) ^ _N[ii](j)) * Mii(k, l);
711 for (int jj = 0; jj != 3; ++jj)
712 if (ii != jj) {
713 Mij(i, j) = (_N[ii](i) ^ _N[jj](j)) + (_N[jj](i) ^ _N[ii](j));
714 P3(i, j, k, l) += ((Vi(ii, jj) * _N[ii](i)) ^ _N[jj](j)) * Mij(k, l);
715 }
716 }
717
718 if (!TANGENT)
719 return P3;
720
721 TLs3(i, j, k, l) = 0.;
722 for (int ii = 0; ii != 3; ++ii)
723 for (int jj = 0; jj != 3; ++jj)
724 Zeta(ii, jj) = stress(i, j) * (_N[ii](i) ^ _N[jj](j));
725
726 for (int ii = 0; ii != 3; ++ii) {
727 Mii(i, j) = (_N[ii](i) ^ _N[ii](j)) + (_N[ii](i) ^ _N[ii](j));
728 TLs3(i, j, k, l) += (0.25 * f(ii) * Zeta(ii, ii) * Mii(i, j)) * Mii(k, l);
729 }
730
731 for (int ii = 0; ii != 3; ++ii)
732 for (int jj = 0; jj != 3; ++jj)
733 if (jj != ii) {
734 Mij(i, j) = (_N[ii](i) ^ _N[jj](j)) + (_N[jj](i) ^ _N[ii](j));
735 Mjj(i, j) = (_N[jj](i) ^ _N[jj](j)) + (_N[jj](i) ^ _N[jj](j));
736 const double k2 = 2. * Ksi(ii, jj);
737 const double zp = k2 * Zeta(ii, jj);
738 const double zp2 = k2 * Zeta(jj, jj);
739 TLs3(i, j, k, l) += (zp * Mij(i, j)) * Mjj(k, l) +
740 (zp * Mjj(i, j)) * Mij(k, l) +
741 (zp2 * Mij(i, j)) * Mij(k, l);
742 for (int kk = 0; kk != 3; ++kk) {
743 if (kk != ii && kk != jj) {
744 Mik(i, j) = (_N[ii](i) ^ _N[kk](j)) + (_N[kk](i) ^ _N[ii](j));
745 Mjk(i, j) = (_N[jj](i) ^ _N[kk](j)) + (_N[kk](i) ^ _N[jj](j));
746 TLs3(i, j, k, l) +=
747 (2. * eta * Zeta(ii, jj) * Mik(i, j)) * Mjk(k, l);
748 }
749 }
750 }
751
752 return P3;
753};
754
755template <typename T>
756inline Tensor2<double, 3, 3>
757getFiniteDiff(const int &kk, const int &ll, CommonData &common_data,
758 Tensor2<T, 3, 3> &F, Tensor2<double, 3, 3> &t_stress) {
760 Tensor2<double, 3, 3> out;
761
762 Tensor2<double, 3, 3> strainPl;
763 Tensor2<double, 3, 3> strainMin;
764 strainPl(i, j) = F(i, j);
765 strainMin(i, j) = F(i, j);
766 const double h = 1e-8;
767 strainPl(kk, ll) += h;
768 strainMin(kk, ll) -= h;
769 Tensor4<PackPtr<double *, 1>, 3, 3, 3, 3> dummy;
770
771 auto calculate_stress = [&](auto &F) {
772 Tensor2<double, 3, 3> Piola;
773 Tensor2_symmetric<double, 3> stress_sym;
774 Tensor2_symmetric<double, 3> C;
775 Tensor2_symmetric<double, 3> strain;
776
777 C(i, j) = F(m, i) ^ F(m, j);
778
780 Tensor2<double, 3, 3> eig_vec;
782
783 auto lnC = get_ln_X_lapack(C, lambda, eig_vec);
784 // auto lnC = get_log_map(C);
785 strain(i, j) = 0.5 * lnC(i, j);
786 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(*common_data.mtD);
787 stress_sym(i, j) = t_D(i, j, k, l) * strain(k, l);
788 auto stress_tmp = to_non_symm(stress_sym);
789
790 // Tensor4<double, 3, 3, 3, 3> D2;
791 // constexpr auto t_kd = Kronecker_Delta<int>();
792 // Tensor4<double, 3, 3, 3, 3> dC_dF;
793 // dC_dF(i, j, k, l) = (t_kd(i, l) * F(k, j)) + (t_kd(j, l) * F(k, i));
794 // auto t_L = get_diff_log_map(C);
795 // t_L(i, j, k, l) *= 0.5;
796 // D2(i, j, k, l) = t_L(i, j, m, n) * dC_dF(m, n, k, l);
797
798 auto D2 = calculate_nominal_moduli_TL(C, F, stress_sym, lambda, eig_vec,
799 dummy, false);
800
801 Piola(k, l) = stress_tmp(i, j) * D2(i, j, k, l);
802
803 return Piola;
804 };
805
806 auto StressPlusH = calculate_stress(strainPl);
807 auto StressMinusH = calculate_stress(strainMin);
808 // out(i, j) = (StressPlusH(i, j) - t_stress(i, j)) / (h);
809
810 out(i, j) = (StressPlusH(i, j) - StressMinusH(i, j)) / (2 * h);
811
812 return out;
813}
814
815template <typename T>
816inline Tensor4<double, 3, 3, 3, 3>
817getDTensorFunction(CommonData &common_data, Tensor2<T, 3, 3> &F,
818 Tensor2<double, 3, 3> &t_stress) {
819 Tensor4<double, 3, 3, 3, 3> my_D;
820 for (int k = 0; k != 3; ++k) {
821 for (int l = 0; l != 3; ++l) {
822 auto d_stress = getFiniteDiff(k, l, common_data, F, t_stress);
823 for (int i = 0; i != 3; ++i) {
824 for (int j = 0; j != 3; ++j) {
825 my_D(i, j, k, l) = d_stress(i, j);
826 }
827 }
828 }
829 }
830 return my_D;
831};
832
833//! [Operators definitions]
834typedef boost::function<Tensor1<double, 3>(const double, const double,
835 const double)>
837
838struct OpStrain : public DomainEleOp {
839 OpStrain(const std::string field_name,
840 boost::shared_ptr<CommonData> common_data_ptr);
841 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
842
843private:
844 boost::shared_ptr<CommonData> commonDataPtr;
845};
846
847struct OpLogStrain : public DomainEleOp {
848 OpLogStrain(const std::string field_name,
849 boost::shared_ptr<CommonData> common_data_ptr);
850 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
851
852private:
853 boost::shared_ptr<CommonData> commonDataPtr;
854};
855
856struct OpStress : public DomainEleOp {
857 OpStress(const std::string field_name,
858 boost::shared_ptr<CommonData> common_data_ptr, boost::shared_ptr<MatrixDouble> mat_D_Ptr);
859 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
860
861private:
862 boost::shared_ptr<CommonData> commonDataPtr;
863 boost::shared_ptr<MatrixDouble> matDPtr;
864};
865
866template <bool TANGENT>
868 OpCalculateLogStressTangent(const std::string field_name,
869 boost::shared_ptr<CommonData> common_data_ptr,
870 boost::shared_ptr<MatrixDouble> mat_D_ptr);
871 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
872
873private:
874 boost::shared_ptr<CommonData> commonDataPtr;
875 boost::shared_ptr<MatrixDouble> matDPtr;
876};
877
880 const std::string field_name,
881 boost::shared_ptr<CommonData> common_data_ptr);
882 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
883
884private:
885 boost::shared_ptr<CommonData> commonDataPtr;
886};
887
888template <bool LOGSTRAIN, bool REACTIONS>
890 OpInternalForceRhs(const std::string field_name,
891 boost::shared_ptr<CommonData> common_data_ptr);
892 MoFEMErrorCode iNtegrate(EntData &data);
893 MoFEMErrorCode aSsemble(EntData &data);
894
895private:
896 boost::shared_ptr<CommonData> commonDataPtr;
897};
898
899template <bool LOGSTRAIN> struct OpPostProcElastic : public DomainEleOp {
900 OpPostProcElastic(const std::string field_name,
901 moab::Interface &post_proc_mesh,
902 std::vector<EntityHandle> &map_gauss_pts,
903 boost::shared_ptr<CommonData> common_data_ptr);
904 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
905
906private:
907 moab::Interface &postProcMesh;
908 std::vector<EntityHandle> &mapGaussPts;
909 boost::shared_ptr<CommonData> commonDataPtr;
910};
911
913 OpSetMaterialBlock(const std::string field_name,
914 boost::shared_ptr<CommonData> common_data_ptr,
915 std::map<EntityHandle, int> &map_block_tets);
916 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
917
918private:
919 boost::shared_ptr<CommonData> commonDataPtr;
920 std::map<EntityHandle, int> mapBlockTets;
921};
922
924
925 OpSetMaterialBlockBoundary(const std::string field_name,
926 boost::shared_ptr<CommonData> common_data_ptr,
927 std::map<EntityHandle, int> &map_block_tets);
928 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
929
930private:
931 boost::shared_ptr<CommonData> commonDataPtr;
932 std::map<EntityHandle, int> mapBlockTets;
933};
934
936 OpSaveReactionForces(MoFEM::Interface &m_field, const std::string field_name,
937 moab::Interface &post_proc_mesh,
938 std::vector<EntityHandle> &map_gauss_pts,
939 boost::shared_ptr<CommonData> common_data_ptr);
940 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
941
942private:
944 moab::Interface &postProcMesh;
945 std::vector<EntityHandle> &mapGaussPts;
946 boost::shared_ptr<CommonData> commonDataPtr;
947};
948
950
951 OpEssentialRHS_U(std::string disp_name,
952 boost::shared_ptr<CommonData> common_data_ptr,
953 VectorDouble &disp, const Range &essential_ents)
954 : BoundaryEleOpAssembly(disp_name, disp_name, BoundaryEleOp::OPROW),
955 essentialEnts(essential_ents), disP(disp),
956 commonDataPtr(common_data_ptr) {}
957
958 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data) {
960 const int nb_row_dofs = row_data.getIndices().size();
961 const int nb_row_base_functions = row_data.getN().size2();
962
963 if (nb_row_dofs) {
964
965 EntityHandle ent = getFEEntityHandle();
966 if (essentialEnts.find(ent) == essentialEnts.end()) {
968 }
969
970 auto t_disp = getFTensor1FromMat<3>(*commonDataPtr->mDispPtr);
971
972 const int nb_integration_pts = getGaussPts().size2();
973 auto t_row_val = row_data.getFTensor0N();
974 auto t_w = getFTensor0IntegrationWeight();
975 const double vol = getMeasure();
976 auto t_coords = getFTensor1CoordsAtGaussPts();
977
978 Tensor1<double, 3> bc_disp(disP(0), disP(1), disP(2));
979 for (int gg = 0; gg != nb_integration_pts; ++gg) {
980 const double a = vol * t_w;
981
982 // // for testing radial displacement
983 // auto theta = atan2(t_coords(1), t_coords(0));
984 // bc_disp(0) = disP(0) * cos(theta);
985 // bc_disp(1) = disP(1) * sin(theta);
986
987 auto t_nf = getFTensor1FromArray<3, 3>(locF);
988 int rr = 0;
989 for (;rr != nb_row_dofs / 3; ++rr) {
990 t_nf(i) += a * t_row_val * (t_disp(i) - bc_disp(i));
991 ++t_row_val;
992 ++t_nf;
993 }
994 for (; rr < nb_row_base_functions; ++rr)
995 ++t_row_val;
996 ++t_disp;
997 ++t_coords;
998 ++t_w;
999 }
1000 }
1002 }
1003
1004private:
1005 boost::shared_ptr<CommonData> commonDataPtr;
1006 VectorDouble &disP;
1008};
1009
1011
1012 OpRotate(std::string disp_name, VectorDouble &omega,
1013 boost::shared_ptr<VectorDouble> c_coords,
1014 boost::shared_ptr<Range> ents_ptr)
1015 : BoundaryEleOpAssembly(disp_name, disp_name, OPROW), entsPtr(ents_ptr),
1016 sourceVec(omega), centerCoords(c_coords) {}
1017
1018 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data) {
1020 const int nb_row_dofs = row_data.getIndices().size();
1021 const int nb_row_base_functions = row_data.getN().size2();
1022
1023 if (nb_row_dofs) {
1024
1025 EntityHandle ent = getFEEntityHandle();
1026 if (entsPtr->find(ent) == entsPtr->end()) {
1028 }
1029
1030 auto t_coords = getFTensor1CoordsAtGaussPts();
1031
1032 const int nb_integration_pts = getGaussPts().size2();
1033 auto t_row_val = row_data.getFTensor0N();
1034 auto t_w = getFTensor0IntegrationWeight();
1035 const double vol = getMeasure();
1037 auto t_omega = getFTensor1FromArray<3, 3>(sourceVec);
1038 auto t_cen_coords = getFTensor1FromArray<3, 3>(*centerCoords);
1040
1041 auto rot = get_rotation_R(t_omega, 1);
1042
1043 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1044 const double a = vol * t_w;
1045 c_dist(i) = t_cen_coords(i) - t_coords(i);
1046 rot_disp(i) = c_dist(i) - rot(j, i) * c_dist(j);
1047
1048 auto t_nf = getFTensor1FromArray<3, 3>(locF);
1049 int rr = 0;
1050 for (; rr != nb_row_dofs / 3; ++rr) {
1051 t_nf(i) -= a * t_row_val * rot_disp(i);
1052 ++t_row_val;
1053 ++t_nf;
1054 }
1055 for (; rr < nb_row_base_functions; ++rr)
1056 ++t_row_val;
1057
1058 ++t_coords;
1059 ++t_w;
1060 }
1061 }
1063 }
1064
1065private:
1066 VectorDouble locRes;
1067 ScalarFun betaCoeff;
1068 boost::shared_ptr<Range> entsPtr;
1069 VectorDouble &sourceVec;
1070 boost::shared_ptr<VectorDouble> centerCoords;
1071};
1072
1073struct OpSetUnSetRowData : public ForcesAndSourcesCore::UserDataOperator {
1074 OpSetUnSetRowData(std::string field_name, bool yes_set,
1075 boost::shared_ptr<CommonData> common_data_ptr)
1076 : ForcesAndSourcesCore::UserDataOperator(
1077 field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
1078 yesSet(yes_set), commonDataPtr(common_data_ptr) {}
1079 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1082 if (!data.getIndices().empty())
1083 if (auto e_ptr = data.getFieldEntities()[0])
1084 if (auto stored_data_ptr =
1085 e_ptr->getSharedStoragePtr<EssentialBcStorage>()) {
1086 if (yesSet) {
1087 commonDataPtr->rowIndices = data.getIndices();
1088 data.getIndices() = stored_data_ptr->entityIndices;
1089 } else {
1090 data.getIndices() = commonDataPtr->rowIndices;
1091 }
1092 }
1093
1095 }
1096
1097public:
1099 boost::shared_ptr<CommonData> commonDataPtr;
1100};
1101
1102struct OpSchurBegin : public DomainEleOp {
1103 OpSchurBegin(const std::string row_field,
1104 boost::shared_ptr<CommonData> common_data_ptr)
1105 : DomainEleOp(row_field, row_field, OPROW),
1106 commonDataPtr(common_data_ptr) {
1107 sYmm = false;
1108 }
1109 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &data);
1110 boost::shared_ptr<CommonData> commonDataPtr;
1111};
1112
1113struct OpSchurEnd : public DomainEleOp {
1114 OpSchurEnd(const std::string row_field,
1115 boost::shared_ptr<CommonData> common_data_ptr,
1116 const double eps = 1e-12)
1117 : DomainEleOp(row_field, row_field, OPROW),
1118 commonDataPtr(common_data_ptr), eps(eps) {
1119 sYmm = false;
1120 }
1121
1122 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
1123
1124private:
1125 const double eps;
1126
1127 MatrixDouble invMat;
1128 MatrixDouble K;
1129 VectorInt iPIV;
1130 VectorDouble lapackWork;
1131
1132 map<std::string, MatrixDouble> invBlockMat;
1133 map<std::string, BlockMatContainer> blockMat;
1134 boost::shared_ptr<CommonData> commonDataPtr;
1135};
1136
1138 OpSchurBeginBoundary(const std::string row_field,
1139 boost::shared_ptr<CommonData> common_data_ptr)
1140 : BoundaryEleOp(row_field, row_field, OPROW),
1141 commonDataPtr(common_data_ptr) {
1142 sYmm = false;
1143 }
1144 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &data);
1145 boost::shared_ptr<CommonData> commonDataPtr;
1146};
1147
1149 OpSchurEndBoundary(const std::string row_field,
1150 boost::shared_ptr<CommonData> common_data_ptr,
1151 const double eps = 1e-12)
1152 : BoundaryEleOp(row_field, row_field, OPROW),
1153 commonDataPtr(common_data_ptr), eps(eps) {
1154 sYmm = false;
1155 }
1156
1157 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
1158
1159private:
1160 const double eps;
1161
1162 MatrixDouble invMat;
1163 MatrixDouble K;
1164 VectorInt iPIV;
1165 VectorDouble lapackWork;
1166
1167 map<std::string, MatrixDouble> invBlockMat;
1168 map<std::string, BlockMatContainer> blockMat;
1169 boost::shared_ptr<CommonData> commonDataPtr;
1170};
1171
1172template <typename TYPE>
1174 OpSchurPreconditionMassContact(const std::string row_field_name,
1175 const std::string col_field_name, ScalarFun beta,
1176 boost::shared_ptr<CommonData> common_data_ptr)
1177 : TYPE(row_field_name, col_field_name, beta),
1178 commonDataPtr(common_data_ptr) {}
1179
1180 MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data,
1181 const bool trans) {
1183 commonDataPtr->massMatTensor = this->locMat;
1185 }
1186
1187private:
1188 boost::shared_ptr<CommonData> commonDataPtr;
1189};
1190
1191// template struct OpSchurPreconditionMassContact<OpMassPrecHTensorBound>;
1192// template struct OpSchurPreconditionMassContact<OpMassPrecHTensorVol>;
1193
1196
1198 OpSchurPreconditionMassTau(const std::string row_field_name,
1199 const std::string col_field_name, ScalarFun beta,
1200 boost::shared_ptr<CommonData> common_data_ptr)
1201 : OpMassPrecScalar(row_field_name, col_field_name, beta),
1202 commonDataPtr(common_data_ptr) {}
1203
1204 MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data,
1205 const bool trans) {
1207 commonDataPtr->massMatScalar = this->locMat;
1209 }
1210
1211private:
1212 boost::shared_ptr<CommonData> commonDataPtr;
1213};
1214
1215
1216/**
1217 * \brief Not used at this stage. Could be used to do some calculations,
1218 * before assembly of local elements.
1219 */
1221
1223 boost::ptr_vector<DataFromMove> &bcData;
1224
1225 boost::shared_ptr<MethodForForceScaling> methodsOpForMove;
1226 boost::shared_ptr<MethodForForceScaling> methodsOpForRollersDisp;
1227
1228 boost::shared_ptr<CommonData> commonDataPtr;
1229
1230 // SmartPetscObj<Mat> SEpEp;
1231 // SmartPetscObj<AO> aoSEpEp;
1232
1234
1236 boost::ptr_vector<DataFromMove> &bc_data,
1237 boost::shared_ptr<CommonData> common_data_ptr)
1238 : mField(m_field), bcData(bc_data), commonDataPtr(common_data_ptr) {}
1239
1240 MoFEMErrorCode preProcess() {
1242 runOnceFlag = true;
1243 CHKERR VecSetOption(ts_F, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
1244
1245 switch (ts_ctx) {
1246 case CTX_TSSETIJACOBIAN: {
1247 snes_ctx = CTX_SNESSETJACOBIAN;
1248 snes_B = ts_B;
1249 // CHKERR MatSetOption(ts_B, MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE);
1250 // CHKERR MatSetOption(ts_B, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
1251
1252 if (commonDataPtr->SEpEp)
1253 CHKERR MatZeroEntries(commonDataPtr->SEpEp);
1254 if (commonDataPtr->STauTau)
1255 CHKERR MatZeroEntries(commonDataPtr->STauTau);
1256
1257 break;
1258 }
1259 case CTX_TSSETIFUNCTION: {
1260 snes_ctx = CTX_SNESSETFUNCTION;
1261 snes_f = ts_F;
1262 break;
1263 }
1264 default:
1265 break;
1266 }
1267
1268 auto time = ts_t;
1269
1270 // forceDataScaled = forceData;
1271 // pressureDataScaled = pressureData;
1272 if (methodsOpForMove)
1273 for (auto &bdata : bcData) {
1274 bdata.scaledValues = bdata.values;
1276 bdata.scaledValues);
1277 }
1278
1279 //FIXME: obsolete, for debugging get rid of this flag
1280 if (time <= (*cache).rollers_stop_time) {
1281
1283 for (auto &roller : (*cache).rollerDataVec) {
1284 roller.BodyDispScaled = roller.rollerDisp;
1286 this, methodsOpForRollersDisp, roller.BodyDispScaled);
1287 }
1288
1289 for (auto &roller : (*cache).rollerDataVec) {
1290 if (roller.methodOpForRollerPosition) {
1291 roller.BodyDispScaled.clear();
1293 this, roller.methodOpForRollerPosition, roller.BodyDispScaled);
1294 }
1295 if (roller.methodOpForRollerDirection) {
1296 roller.BodyDirectionScaled.clear();
1298 this, roller.methodOpForRollerDirection,
1299 roller.BodyDirectionScaled);
1300 }
1301 }
1302 }
1303
1304 // if (snes_ctx == CTX_SNESNONE) {
1305 // if (!dofsIndices.empty()) {
1306 // CHKERR VecSetValues(snes_x, dofsIndices.size(),
1307 // &*dofsIndices.begin(),
1308 // &*dofsValues.begin(), INSERT_VALUES);
1309 // }
1310 // CHKERR VecAssemblyBegin(snes_x);
1311 // CHKERR VecAssemblyEnd(snes_x);
1312 // }
1313
1315 }
1316
1317 MoFEMErrorCode postProcess() {
1319 // if (ts_ctx == CTX_TSSETIJACOBIAN)
1320 // CHKERR MatView(ts_B, PETSC_VIEWER_STDOUT_SELF);
1321 auto assemble = [](Mat a) {
1323 if (a) {
1324 CHKERR MatAssemblyBegin(a, MAT_FINAL_ASSEMBLY);
1325 CHKERR MatAssemblyEnd(a, MAT_FINAL_ASSEMBLY);
1326 }
1328 };
1329
1330 if (ts_ctx == CTX_TSSETIJACOBIAN) {
1331
1332 if (commonDataPtr->SEpEp)
1333 CHKERR assemble(commonDataPtr->SEpEp);
1334 if (commonDataPtr->STauTau)
1335 CHKERR assemble(commonDataPtr->STauTau);
1336
1337 // if (true)
1338 // CHKERR MatView(ts_B, PETSC_VIEWER_DRAW_WORLD);
1339 // if (commonDataPtr->SEpEp)
1340 // CHKERR MatView(commonDataPtr->SEpEp, PETSC_VIEWER_DRAW_WORLD);
1341 // if (commonDataPtr->STauTau)
1342 // CHKERR MatView(commonDataPtr->STauTau, PETSC_VIEWER_DRAW_WORLD);
1343
1344 }
1345
1347 }
1348 };
1349}; // namespace OpElasticTools
static Number< 2 > N2
static Number< 1 > N1
static Number< 0 > N0
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
static const double eps
Kronecker Delta class.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
constexpr double omega
Save field DOFS on vertices/tags.
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
@ F
constexpr auto t_kd
double eta
FTensor::Index< 'i', SPACE_DIM > i
static double lambda
const double v
phase velocity of light in medium (cm/ns)
static __CLPK_integer lapack_dsyev(char jobz, char uplo, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *w, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition: lapack_wrap.h:261
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1884
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const double T
FTensor::Ddg< double, 3, 3 > getDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, const int nb)
Get the Diff Mat object.
FTensor::Ddg< double, 3, 3 > getDiffDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, Fun< double > dd_f, FTensor::Tensor2< double, 3, 3 > &t_S, const int nb)
MoFEMErrorCode get_eigen_val_and_proj_lapack(const Tensor2_symmetric< T, 3 > &X, Tensor1< T, 3 > &lambda, Tensor2< T, 3, 3 > &eig_vec)
auto get_ln_X(const Tensor2_symmetric< T, 3 > &X)
bool is_ddg(Tensor4< T, 3, 3, 3, 3 > tensor)
bool is_diff(const double &a, const double &b)
Ddg< double, 3, 3 > tensor4_to_ddg(Tensor4< T, 3, 3, 3, 3 > &t)
Tensor2< double, 3, 3 > getFiniteDiff(const int &kk, const int &ll, CommonData &common_data, Tensor2< T, 3, 3 > &F, Tensor2< double, 3, 3 > &t_stress)
auto calculate_lagrangian_moduli_TL(Tensor2_symmetric< TP, 3 > &stress, Tensor1< T, 3 > &lambda, Tensor2< T, 3, 3 > &eig_vec, Ddg< T3, 3, 3 > &TLs3)
boost::function< Tensor1< double, 3 >(const double, const double, const double)> VectorFun
[Operators definitions]
auto get_ln_X_lapack(const Tensor2_symmetric< T1, 3 > &X, Tensor1< T2, 3 > &lambda, Tensor2< T3, 3, 3 > &eig_vec)
auto get_log_map(const Tensor2_symmetric< T, 3 > &X)
auto get_tensor4_from_ddg(Ddg< T, 3, 3 > &ddg)
Tensor4< double, 3, 3, 3, 3 > getDTensorFunction(CommonData &common_data, Tensor2< T, 3, 3 > &F, Tensor2< double, 3, 3 > &t_stress)
auto get_diff_log_map(Tensor2_symmetric< T, 3 > &X)
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMass< 3, 9 > OpMassPrecHTensorBound
auto to_non_symm(T &symm)
auto get_rotation_R(Tensor1< T, 3 > t1_omega, double tt)
bool is_eq(const double &a, const double &b)
auto to_symm(T &non_symm)
auto exp_map(const Tensor2_symmetric< T, 3 > &X)
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMass< 3, 9 > OpMassPrecHTensorVol
auto calculate_nominal_moduli_TL(Tensor2_symmetric< T, 3 > &X, Tensor2< TP, 3, 3 > &F, Tensor2_symmetric< TP, 3 > &stress, Tensor1< T, 3 > &lambda, Tensor2< T, 3, 3 > &eig_vec, Tensor4< T3, 3, 3, 3, 3 > &TLs3, bool compute_tangent)
auto test_projection_lukasz(Tensor2_symmetric< T, 3 > &X, Tensor2_symmetric< TP, 3 > &stress, Tensor1< T, 3 > &lambda, Tensor2< T, 3, 3 > &eig_vec, Ddg< T3, 3, 3 > &TLs3)
double h
constexpr double t
plate stiffness
Definition: plate.cpp:59
constexpr auto field_name
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
bool sYmm
If true assume that matrix is symmetric structure.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
const VectorFieldEntities & getFieldEntities() const
get field entities
const VectorInt & getIndices() const
Get global indices of dofs on entity.
@ OPROW
operator doWork function is executed on FE rows
array< double, 2 > stateArrayArea
boost::shared_ptr< MatrixDouble > mPiolaStressPtr
SmartPetscObj< Mat > STauTau
SmartPetscObj< AO > aoSTauTau
SmartPetscObj< Vec > reactionsVec
boost::shared_ptr< VectorDouble > detF
boost::shared_ptr< BlockMatContainer > blockMatContainerPtr
SmartPetscObj< Mat > SEpEp
boost::shared_ptr< MatrixDouble > mStrainPtr
boost::shared_ptr< MatrixDouble > matEigenVal
boost::shared_ptr< MatrixDouble > mStressPtr
SmartPetscObj< AO > aoSEpEp
boost::shared_ptr< MatrixDouble > mtD_Deviator
boost::shared_ptr< MatrixDouble > mGradPtr
boost::shared_ptr< MatrixDouble > mtD_Axiator
boost::shared_ptr< MatrixDouble > materialTangent
boost::shared_ptr< MatrixDouble > mDeformationGradient
boost::shared_ptr< MatrixDouble > mDispPtr
boost::shared_ptr< MatrixDouble > matC
boost::shared_ptr< MatrixDouble > matEigenVector
boost::shared_ptr< MatrixDouble > dE_dF_mat
boost::shared_ptr< MatrixDouble > mtD
Not used at this stage. Could be used to do some calculations, before assembly of local elements.
MoFEMErrorCode preProcess()
boost::shared_ptr< MethodForForceScaling > methodsOpForMove
boost::ptr_vector< DataFromMove > & bcData
MoFEMErrorCode postProcess()
bool runOnceFlag
FePrePostProcess(MoFEM::Interface &m_field, boost::ptr_vector< DataFromMove > &bc_data, boost::shared_ptr< CommonData > common_data_ptr)
MoFEM::Interface & mField
boost::shared_ptr< MethodForForceScaling > methodsOpForRollersDisp
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< MatrixDouble > matDPtr
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< CommonData > commonDataPtr
OpEssentialRHS_U(std::string disp_name, boost::shared_ptr< CommonData > common_data_ptr, VectorDouble &disp, const Range &essential_ents)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data)
MoFEMErrorCode iNtegrate(EntData &data)
[Internal force]
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode aSsemble(EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Calculate strain]
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Postprocessing]
std::vector< EntityHandle > & mapGaussPts
moab::Interface & postProcMesh
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data)
OpRotate(std::string disp_name, VectorDouble &omega, boost::shared_ptr< VectorDouble > c_coords, boost::shared_ptr< Range > ents_ptr)
boost::shared_ptr< Range > entsPtr
boost::shared_ptr< VectorDouble > centerCoords
std::vector< EntityHandle > & mapGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< CommonData > commonDataPtr
OpSchurBeginBoundary(const std::string row_field, boost::shared_ptr< CommonData > common_data_ptr)
OpSchurBegin(const std::string row_field, boost::shared_ptr< CommonData > common_data_ptr)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
map< std::string, MatrixDouble > invBlockMat
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
map< std::string, BlockMatContainer > blockMat
OpSchurEndBoundary(const std::string row_field, boost::shared_ptr< CommonData > common_data_ptr, const double eps=1e-12)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpSchurEnd(const std::string row_field, boost::shared_ptr< CommonData > common_data_ptr, const double eps=1e-12)
boost::shared_ptr< CommonData > commonDataPtr
map< std::string, MatrixDouble > invBlockMat
map< std::string, BlockMatContainer > blockMat
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data, const bool trans)
OpSchurPreconditionMassContact(const std::string row_field_name, const std::string col_field_name, ScalarFun beta, boost::shared_ptr< CommonData > common_data_ptr)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data, const bool trans)
OpSchurPreconditionMassTau(const std::string row_field_name, const std::string col_field_name, ScalarFun beta, boost::shared_ptr< CommonData > common_data_ptr)
std::map< EntityHandle, int > mapBlockTets
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< CommonData > commonDataPtr
std::map< EntityHandle, int > mapBlockTets
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpSetUnSetRowData(std::string field_name, bool yes_set, boost::shared_ptr< CommonData > common_data_ptr)
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Calculate strain]
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Calculate stress]
boost::shared_ptr< MatrixDouble > matDPtr