v0.13.2
Loading...
Searching...
No Matches
ElasticOperators.cpp
Go to the documentation of this file.
1/** \file ElasticOperators.cpp
2 */
3
4/* This file is part of MoFEM.
5 * MoFEM is free software: you can redistribute it and/or modify it under
6 * the terms of the GNU Lesser General Public License as published by the
7 * Free Software Foundation, either version 3 of the License, or (at your
8 * option) any later version.
9 *
10 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
11 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
13 * License for more details.
14 *
15 * You should have received a copy of the GNU Lesser General Public
16 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
17
18#include <MoFEM.hpp>
19using namespace MoFEM;
20using namespace FTensor;
21
24#include <RigidBodies.hpp>
25#include <BlockMatData.hpp>
27#include <BasicFeTools.hpp>
28#include <MatrixFunction.hpp>
29
30#include <ElasticOperators.hpp>
31#include <chrono> //for debugging
32
33struct MeasureTime {
34 int time;
35 chrono::high_resolution_clock::time_point start;
36 chrono::high_resolution_clock::time_point stop;
37 MeasureTime() { start = chrono::high_resolution_clock::now(); }
39 stop = chrono::high_resolution_clock::now();
40 auto duration = chrono::duration_cast<chrono::microseconds>(stop - start);
41 cout << "Time taken by function: " << duration.count() << " ms." << endl;
42 }
43};
44
45namespace OpElasticTools {
46
48 boost::shared_ptr<CommonData> common_data_ptr)
50 commonDataPtr(common_data_ptr) {
51 // Operator is only executed for vertices
52 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
53}
54
55//! [Calculate strain]
58 const size_t nb_gauss_pts = commonDataPtr->mGradPtr->size2();
59 commonDataPtr->mStrainPtr->resize(6, nb_gauss_pts);
60 auto t_grad = getFTensor2FromMat<3, 3>(*(commonDataPtr->mGradPtr));
61 auto t_strain = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStrainPtr));
62 commonDataPtr->detF->resize(nb_gauss_pts, false);
63 auto t_detF = getFTensor0FromVec(*(commonDataPtr->detF));
64
65 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
66 t_strain(i, j) = (t_grad(i, j) || t_grad(j, i)) / 2;
67 t_detF = 1;
68 ++t_detF;
69 ++t_grad;
70 ++t_strain;
71 }
72
74}
75
77 boost::shared_ptr<CommonData> common_data_ptr)
79 commonDataPtr(common_data_ptr) {
80 // Operator is only executed for vertices
81 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
82}
83
84//! [Calculate strain]
87 const size_t nb_gauss_pts = commonDataPtr->mGradPtr->size2();
88 commonDataPtr->mStrainPtr->resize(6, nb_gauss_pts);
89 commonDataPtr->mDeformationGradient->resize(9, nb_gauss_pts, false);
90 commonDataPtr->matC->resize(6, nb_gauss_pts, false);
91 commonDataPtr->matEigenVal->resize(3, nb_gauss_pts, false);
92 commonDataPtr->matEigenVector->resize(9, nb_gauss_pts, false);
93 commonDataPtr->detF->resize(nb_gauss_pts, false);
94
95 // double dt = 0.;
96 // auto get_dt = [&]() {
97 // double dt;
98 // CHKERR TSGetTime(getFEMethod()->ts, &dt);
99 // // CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
100 // dt = dt;
101 // return dt;
102 // };
103 // if ((getDataCtx() & PetscData::CtxSetTime).any())
104 // dt = get_dt();
105
106 auto t_grad = getFTensor2FromMat<3, 3>(*(commonDataPtr->mGradPtr));
107 auto F = getFTensor2FromMat<3, 3>(*(commonDataPtr->mDeformationGradient));
108 auto C = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->matC));
109 auto t_eigen = getFTensor1FromMat<3>(*(commonDataPtr->matEigenVal));
110 auto t_eigen_vec = getFTensor2FromMat<3, 3>(*(commonDataPtr->matEigenVector));
111 auto t_detF = getFTensor0FromVec(*(commonDataPtr->detF));
112 auto t_strain = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStrainPtr));
113
114 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
115 // auto R = get_rotation_R((*cache).omega, dt);
116 // F(i, j) = (t_grad(i, k) + kronecker_delta(i, k)) * R(k, j);
117 F(i, j) = t_grad(i, j) + kronecker_delta(i, j);
118
119 double det;
121 t_detF = det;
122
123 // test
124 // if (getTStime() > 0.2) {
125 // F(k, i) = 0.;
126 // F(0, 0) = 1;
127 // F(1, 1) = 1;
128 // F(2, 2) = 3;
129 // F(0, 1) = 2;
130 // F(1, 0) = 2;
131 // }
132
133 C(i, j) = F(k, i) ^ F(k, j);
134
135 CHKERR get_eigen_val_and_proj_lapack(C, t_eigen, t_eigen_vec);
136 auto lnC = get_ln_X_lapack(C, t_eigen, t_eigen_vec);
137 t_strain(i, j) = 0.5 * lnC(i, j);
138
139 ++t_eigen;
140 ++t_eigen_vec;
141 ++C;
142 ++t_detF;
143 ++F;
144 ++t_grad;
145 ++t_strain;
146 }
147
149}
151 boost::shared_ptr<CommonData> common_data_ptr, boost::shared_ptr<MatrixDouble> mat_D_Ptr)
153 commonDataPtr(common_data_ptr), matDPtr(mat_D_Ptr) {
154 // Operator is only executed for vertices
155 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
156}
157
158//! [Calculate stress]
161 const size_t nb_gauss_pts = commonDataPtr->mStrainPtr->size2();
162 commonDataPtr->mStressPtr->resize(6, nb_gauss_pts);
163 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(*matDPtr);
164 auto t_strain = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStrainPtr));
165 auto t_stress = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStressPtr));
166
167 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
168 t_stress(i, j) = t_D(i, j, k, l) * t_strain(k, l);
169 ++t_strain;
170 ++t_stress;
171 }
172
174}
175template <bool TANGENT>
177 const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr,
178 boost::shared_ptr<MatrixDouble> mat_D_ptr)
180 commonDataPtr(common_data_ptr), matDPtr(mat_D_ptr) {
181 // Operator is only executed for vertices
182 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
183}
184
185template <bool TANGENT>
187 EntityType type,
188 EntData &data) {
190
191 const size_t nb_gauss_pts = data.getN().size1();
192 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(*matDPtr);
193 auto t_stress = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStressPtr));
194 commonDataPtr->mPiolaStressPtr->resize(9, nb_gauss_pts, false);
195 auto t_piola = getFTensor2FromMat<3, 3>(*(commonDataPtr->mPiolaStressPtr));
196 auto F = getFTensor2FromMat<3, 3>(*(commonDataPtr->mDeformationGradient));
197 auto C = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->matC));
198 auto t_eigen_vec = getFTensor2FromMat<3, 3>(*(commonDataPtr->matEigenVector));
199 auto t_eigen = getFTensor1FromMat<3>(*(commonDataPtr->matEigenVal));
200
201 MatrixDouble &dP_dF = *(commonDataPtr->materialTangent);
202 dP_dF.resize(81, nb_gauss_pts, false);
203 auto t_dP_dF = getFTensor4FromMat<3, 3, 3, 3>(dP_dF);
204
205 MatrixDouble &dE_dF = *(commonDataPtr->dE_dF_mat);
206 dE_dF.resize(81, nb_gauss_pts, false);
207 auto t_dE_dF = getFTensor4FromMat<3, 3, 3, 3>(dE_dF);
208
209 constexpr auto t_kd = Kronecker_Delta<int>();
211 // Ddg<double, 3, 3> TL_test;
213
214 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
216 dC_dF(i, j, k, l) = (t_kd(i, l) * F(k, j)) + (t_kd(j, l) * F(k, i));
217 auto dlogC_dC = calculate_lagrangian_moduli_TL<TANGENT>(t_stress, t_eigen,
218 t_eigen_vec, TL);
219
220 t_dE_dF(i, j, k, l) = dlogC_dC(i, j, m, n) * dC_dF(m, n, k, l);
221 t_dE_dF(i, j, k, l) *= 0.5;
222
223 if (TANGENT) {
224 S(k, l) = t_stress(i, j) * dlogC_dC(i, j, k, l);
225 t_piola(i, l) = F(i, k) * S(k, l);
226
227 FTensor::Ddg<double, 3, 3> P_D_P_plus_TL;
228 P_D_P_plus_TL(i, j, k, l) =
229 TL(i, j, k, l) +
230 (dlogC_dC(i, j, o, p) * t_D(o, p, m, n)) * dlogC_dC(m, n, k, l);
231 P_D_P_plus_TL(i, j, k, l) *= 0.5;
232
233 t_dP_dF(i, j, m, n) = t_kd(i, m) * (t_kd(k, n) * S(k, j));
234 t_dP_dF(i, j, m, n) +=
235 F(i, k) * (P_D_P_plus_TL(k, j, o, p) * dC_dF(o, p, m, n));
236
237 // auto Piola = to_non_symm(t_stress);
238 // auto D1_from_fddsa321p = getDTensorFunction(*commonDataPtr, F, Piola);
239 // auto test_lukasz =
240 // test_projection_lukasz(C, t_stress, t_eigen, t_eigen_vec, TL_test);
241 // test_lukasz(i, j, k, l) *= 2;
242 // TL_test(i, j, k, l) *= 4;
243 // string wait;
244 } else {
245 // auto dlogC_dC = calculate_lagrangian_moduli_TL<false>(
246 // t_stress, t_eigen, t_eigen_vec, TL);
247 S(k, l) = t_stress(i, j) * dlogC_dC(i, j, k, l);
248 t_piola(i, l) = F(i, k) * S(k, l);
249 }
250 if (false && TANGENT) {
251 if (getTStime() > 0.2) {
252 {
253 Ddg<double, 3, 3> TL3_lukasz;
254 cout << "NEW IMPLEMENTATION" << endl;
255 MeasureTime time1;
256 for (int ii = 0; ii != 1000; ii++)
257 auto P_proj_lukasz = test_projection_lukasz(
258 C, t_stress, t_eigen, t_eigen_vec, TL3_lukasz);
259 }
260 Ddg<double, 3, 3> dummy33;
261 {
262 cout << "NEW IMPLEMENTATION" << endl;
263 MeasureTime time1;
264 for (int ii = 0; ii != 1000; ii++)
265 auto P_proj_ = calculate_lagrangian_moduli_TL<true>(
266 t_stress, t_eigen, t_eigen_vec, dummy33);
267 }
268 {
269 cout << "(OLD) FINITE DIFFERENCE" << endl;
270 auto Piola = to_non_symm(t_stress);
271 MeasureTime time1;
272 for (int ii = 0; ii != 1000; ii++)
273 auto D1_from_fddsap = getDTensorFunction(*commonDataPtr, F, Piola);
274 }
275 {
276 cout << "EIGENVALUES (LAPACK)" << endl;
277 MeasureTime time1;
278 for (int ii = 0; ii != 1000; ii++)
279 CHKERR get_eigen_val_and_proj_lapack(C, t_eigen, t_eigen_vec);
280 }
281 }
282
283 // double sum = 0;
284 // double sum2 = 0;
285 // double sum3 = 0;
286 // double sum0 = 0;
287 // for (int i = 0; i != 3; ++i)
288 // for (int j = 0; j != 3; ++j)
289 // for (int k = 0; k != 3; ++k)
290 // for (int l = 0; l != 3; ++l) {
291 // // sum0 += abs(P3_prop_lag(i, j, k, l) - dlogC_dC_test(i, j, k,
292 // l));
293 // // sum += abs(P_Proj_lukasz(i, j, k, l) - dlogC_dC_test(i, j,
294 // k, l)); sum2 += abs(dP_dF(i, j, k, l) - D1(i, j, k, l)); sum3
295 // +=
296 // abs(myDFin_lukasz(i, j, k, l) - D1_from_fd_LAG(i, j, k,
297 // l));
298 // // if (abs(myDFinFinFin(i, j, k, l) - D1(i, j, k, l)) > 1e-6)
299 // // cout << i << j << k << l << " "
300 // // << myDFinFinFin(i, j, k, l) / D1(i, j, k, l) << endl;
301 // }
302 if (getTStime() > 0.2)
303 string wait;
304 }
305
306 ++t_piola;
307 ++t_eigen;
308 ++t_eigen_vec;
309 ++t_dP_dF;
310 ++t_dE_dF;
311 ++t_stress;
312 ++F;
313 ++C;
314 }
315
317}
318
320 const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
322 commonDataPtr(common_data_ptr) {
323 // Operator is only executed for vertices
324 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
325}
327 EntityType type,
328 EntData &data) {
330
331 const size_t nb_gauss_pts = data.getN().size1();
332 auto t_stress = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStressPtr));
333 commonDataPtr->mPiolaStressPtr->resize(9, nb_gauss_pts, false);
334 auto t_piola = getFTensor2FromMat<3, 3>(*(commonDataPtr->mPiolaStressPtr));
335 auto F = getFTensor2FromMat<3, 3>(*(commonDataPtr->mDeformationGradient));
336 auto t_detF = getFTensor0FromVec(*(commonDataPtr->detF));
337
338 MatrixDouble &dP_dF = *(commonDataPtr->materialTangent);
339 dP_dF.resize(81, nb_gauss_pts, false);
340 auto t_dP_dF = getFTensor4FromMat<3, 3, 3, 3>(dP_dF);
341
342 MatrixDouble &dE_dF = *(commonDataPtr->dE_dF_mat);
343 dE_dF.resize(81, nb_gauss_pts, false);
344 auto t_dE_dF = getFTensor4FromMat<3, 3, 3, 3>(dE_dF);
345 const bool is_lhs = (int)getTSCtx() == 3;
346 constexpr auto t_kd = Kronecker_Delta<int>();
347
348 const double lambda = LAMBDA((*cache).young_modulus, (*cache).poisson);
349 const double mu = MU((*cache).young_modulus, (*cache).poisson);
350
351 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
352
354 CHKERR invertTensor3by3(F, t_detF, invF);
355 const double log_det = log(t_detF);
356
357 if (is_lhs) {
358 // compute tangent
359 const double coef = mu - lambda * log_det;
360 t_dP_dF(i, J, k, L) =
361 invF(J, i) * (invF(L, k) * lambda) + invF(L, i) * (invF(J, k) * coef);
362
363 t_dP_dF(0, 0, 0, 0) += mu;
364 t_dP_dF(0, 1, 0, 1) += mu;
365 t_dP_dF(0, 2, 0, 2) += mu;
366 t_dP_dF(1, 0, 1, 0) += mu;
367 t_dP_dF(1, 1, 1, 1) += mu;
368 t_dP_dF(1, 2, 1, 2) += mu;
369 t_dP_dF(2, 0, 2, 0) += mu;
370 t_dP_dF(2, 1, 2, 1) += mu;
371 t_dP_dF(2, 2, 2, 2) += mu;
372
373 } else {
374
375 // compute only stress
376 const double var = lambda * log_det - mu;
377 t_piola(i, J) = mu * F(i, J) + var * invF(J, i);
378 }
379
380 ++t_detF;
381 ++t_piola;
382 ++t_dP_dF;
383 ++t_stress;
384 ++F;
385 }
386
388}
389
390//! [Calculate stress]
391template <bool LOGSTRAIN, bool REACTIONS>
393 const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
395 commonDataPtr(common_data_ptr) {}
396
397//! [Internal force]
398template <bool LOGSTRAIN, bool REACTIONS>
402
403 const size_t nb_dofs = data.getIndices().size();
404 if (nb_dofs) {
405
406 const size_t nb_base_functions = data.getN().size2();
407 if (3 * nb_base_functions < nb_dofs)
408 SETERRQ(
409 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
410 "Number of DOFs is larger than number of base functions on entity");
411
412 const size_t nb_gauss_pts = data.getN().size1();
413
414 auto t_w = getFTensor0IntegrationWeight();
415 auto t_stress =
416 getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStressPtr));
417 auto t_piola = getFTensor2FromMat<3, 3>(*(commonDataPtr->mPiolaStressPtr));
418 auto t_diff_base = data.getFTensor1DiffN<3>();
419
420 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
421
422 double alpha = getMeasure() * t_w;
423 auto t_nf = getFTensor1FromArray<3, 3>(locF);
424 auto out_stress = to_non_symm(t_stress);
425
426 if (LOGSTRAIN) {
427 out_stress(i, j) = t_piola(i, j);
428 ++t_piola;
429 }
430
431 size_t bb = 0;
432 for (; bb != nb_dofs / 3; ++bb) {
433 t_nf(i) += alpha * t_diff_base(j) * out_stress(i, j);
434 ++t_diff_base;
435 ++t_nf;
436 }
437 for (; bb < nb_base_functions; ++bb)
438 ++t_diff_base;
439
440 ++t_stress;
441 ++t_w;
442 }
443 }
445}
446
447template <bool LOGSTRAIN, bool REACTIONS>
451 if (REACTIONS) {
452 auto react_ents = commonDataPtr->reactionEnts;
453 EntityHandle ent = getFEEntityHandle();
454
455 auto &vec = commonDataPtr->reactionsVec;
456 double *react_ptr;
457 CHKERR VecGetArray(vec, &react_ptr);
458
459 int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
460 int ii = 0;
461 for (auto &dof : data.getFieldDofs()) {
462 auto ent = dof->getEnt();
463 double tag_val[rank];
464 auto idx = dof->getDofCoeffIdx();
465 auto &m_field = getPtrFE()->mField;
466 CHKERR m_field.get_moab().tag_get_data(commonDataPtr->reactionTag, &ent,
467 1, tag_val);
468 if (dof->getEntType() == MBVERTEX) {
469 const double react = this->locF(ii) / (*cache).young_modulus_inv;
470 tag_val[idx] += react;
471 if (react_ents.find(ent) != react_ents.end())
472 react_ptr[idx] += react;
473 }
474 ii++;
475 CHKERR m_field.get_moab().tag_set_data(commonDataPtr->reactionTag, &ent,
476 1, &tag_val);
477 }
478 CHKERR VecRestoreArray(vec, &react_ptr);
479 } else
481 this->getKSPf(), data, &*this->locF.data().begin(), ADD_VALUES);
482
484}
485
487 const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr,
488 std::map<EntityHandle, int> &map_block_tets)
490 commonDataPtr(common_data_ptr), mapBlockTets(map_block_tets) {
491 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
492}
493
495 EntData &data) {
497 if (mapBlockTets.empty())
499
500 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
501 const auto id = mapBlockTets.at(fe_ent);
502 cache = &mat_blocks.at(id);
503 *commonDataPtr->mtD = (*cache).mtD;
504
505 // TODO:
506 // FIXME: include params for plasticity
508}
509
511 const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr,
512 std::map<EntityHandle, int> &map_block_tets)
514 commonDataPtr(common_data_ptr), mapBlockTets(map_block_tets) {
515 std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
516 doEntities[MBTRI] = doEntities[MBQUAD] = true;
517}
518
520 EntData &data) {
522 if (mapBlockTets.empty())
524 std::vector<EntityHandle> adj_ents;
525 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
526 auto &moab = getPtrFE()->mField.get_moab();
527 CHKERR moab.get_adjacencies(&fe_ent, 1, 3, false, adj_ents);
528
529 const auto id = mapBlockTets.at(adj_ents[0]);
530 cache = &mat_blocks.at(id);
531 *commonDataPtr->mtD = (*cache).mtD;
532
534}
535
537 MoFEM::Interface &m_field, const std::string field_name,
538 moab::Interface &post_proc_mesh, std::vector<EntityHandle> &map_gauss_pts,
539 boost::shared_ptr<CommonData> common_data_ptr)
540 : DomainEleOp(field_name, field_name, DomainEleOp::OPROW), mField(m_field),
541 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
542 commonDataPtr(common_data_ptr) {
543 // Operator is only executed for vertices
544 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
545}
546
548 EntData &data) {
550
551 auto get_tag = [&](const std::string name, size_t size) {
552 std::array<double, 9> def;
553 std::fill(def.begin(), def.end(), 0);
554 Tag th;
555 CHKERR postProcMesh.tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE, th,
556 MB_TAG_CREAT | MB_TAG_SPARSE,
557 def.data());
558 return th;
559 };
560
561 auto th_react = get_tag("REACTIONS", 3);
562 size_t nb_gauss_pts = data.getN().size1();
563 VectorDouble vAlues(data.getFieldData().size());
564
565 const size_t rank = data.getFieldDofs()[0]->getNbOfCoeffs();
566
567 int dd = 0;
568 for (auto &dof : data.getFieldDofs()) {
569 auto ent = dof->getEnt();
570 // ents_vec[dd] = ent;
571 const size_t idx = dof->getDofCoeffIdx();
572
573 double tag_val[rank];
574 CHKERR mField.get_moab().tag_get_data(commonDataPtr->reactionTag, &ent, 1,
575 tag_val);
576 vAlues(dd++) = tag_val[idx];
577 }
578
579 const void *tags_ptr[mapGaussPts.size()];
580 CHKERR postProcMesh.tag_get_by_ptr(th_react, &mapGaussPts[0],
581 mapGaussPts.size(), tags_ptr);
582 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
583 for (int rr = 0; rr != rank; rr++) {
584 // its ugly, but tested in PostProcOnRefMesh
585 const double val = cblas_ddot((vAlues.size() / rank), &(data.getN(gg)[0]),
586 1, &vAlues(rr), rank);
587 ((double *)tags_ptr[gg])[rr] += val;
588 }
589 }
590
592}
593template <bool LOGSTRAIN>
595 const std::string field_name, moab::Interface &post_proc_mesh,
596 std::vector<EntityHandle> &map_gauss_pts,
597 boost::shared_ptr<CommonData> common_data_ptr)
599 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
600 commonDataPtr(common_data_ptr) {
601 // Operator is only executed for vertices
602 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
603}
604
605//! [Postprocessing]
606template <bool LOGSTRAIN>
608 EntData &data) {
610
611 auto get_tag = [&](const std::string name, size_t size = 9) {
612 std::array<double, 9> def;
613 std::fill(def.begin(), def.end(), 0);
614 Tag th;
615 CHKERR postProcMesh.tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE, th,
616 MB_TAG_CREAT | MB_TAG_SPARSE,
617 def.data());
618 return th;
619 };
620
621 MatrixDouble3by3 mat(3, 3);
622
623 auto set_matrix = [&](auto &t) -> MatrixDouble3by3 & {
624 mat.clear();
625 for (size_t r = 0; r != 3; ++r)
626 for (size_t c = 0; c != 3; ++c)
627 mat(r, c) = t(r, c);
628 return mat;
629 };
630
631 auto set_scalar = [&](auto t) -> MatrixDouble3by3 & {
632 mat.clear();
633 mat(0, 0) = t;
634 return mat;
635 };
636
637 auto set_tag = [&](auto th, auto gg, MatrixDouble3by3 &mat) {
638 return postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
639 &*mat.data().begin());
640 };
641
642 string strain_tag_name = "STRAIN";
643 string stress_tag_name = "CAUCHY_STRESS";
644 if (LOGSTRAIN) {
645 strain_tag_name = "LOG_STRAIN";
646 // stress_tag_name = "PIOLA_KIRCHHOFF_STRESS";
647 }
648
649 auto th_strain = get_tag(strain_tag_name);
651 auto th_grad = get_tag("GRAD");
652 auto th_stress = get_tag(stress_tag_name);
653 auto th_mises = get_tag("MISES", 1);
654
655 size_t nb_gauss_pts = data.getN().size1();
656 auto t_grad = getFTensor2FromMat<3, 3>(*(commonDataPtr->mGradPtr));
657 auto t_strain = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStrainPtr));
658 auto t_stress = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStressPtr));
659 auto t_piola = getFTensor2FromMat<3, 3>(*(commonDataPtr->mPiolaStressPtr));
660
661 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
662
663 t_stress(i, j) /= (*cache).young_modulus_inv;
664
665 if (LOGSTRAIN) {
666 // auto o_stress = to_non_symm(t_stress);
667 t_piola(i, j) /= (*cache).young_modulus_inv;
668 F(i, j) = t_grad(i, j) + kronecker_delta(i, j);
669 double detF;
671 // cauchy
672 t_stress(i, j) = (1. / detF) * (t_piola(i, k) ^ F(j, k));
673 ++t_piola;
674 }
675
676 const double mises =
677 sqrt(0.5 * (pow(t_stress(0, 0) - t_stress(1, 1), 2) +
678 pow(t_stress(1, 1) - t_stress(2, 2), 2) +
679 pow(t_stress(2, 2) - t_stress(0, 0), 2) +
680 6 * (pow(t_stress(0, 1), 2) + pow(t_stress(1, 2), 2) +
681 pow(t_stress(2, 0), 2))));
682
683 CHKERR set_tag(th_mises, gg, set_scalar(mises));
684 // FIXME: add option for more postprocessing
685 // CHKERR set_tag(th_grad, gg, set_matrix(t_grad));
686 if (!commonDataPtr->isNeoHook)
687 CHKERR set_tag(th_strain, gg, set_matrix(t_strain));
688 CHKERR set_tag(th_stress, gg, set_matrix(t_stress));
689
690 ++t_grad;
691 ++t_strain;
692 ++t_stress;
693 }
694
696}
697
699 EntData &data) {
701 if (row_type != MBTRI && row_type != MBQUAD)
703
704 // MoFEMFunctionReturnHot(0); //FIXME:
705 if (!commonDataPtr->blockMatContainerPtr)
707 auto &bmc = *commonDataPtr->blockMatContainerPtr;
708 for (auto &bit : bmc)
709 bit.unSetAtElement();
711}
712
714 EntData &data) {
716 if (row_type != MBTRI && row_type != MBQUAD)
718
719 // MoFEMFunctionReturnHot(0); //FIXME:
720
721 if (commonDataPtr->blockMatContainerPtr) {
722
723 auto SEpEp = commonDataPtr->SEpEp;
724 auto aoSEpEp = commonDataPtr->aoSEpEp;
725 auto STauTau = commonDataPtr->STauTau;
726 auto aoSTauTau = commonDataPtr->aoSTauTau;
727
728 if (SEpEp) {
729
730 auto find_block = [&](BlockMatContainer &add_bmc, auto &row_name,
731 auto &col_name, auto row_side, auto col_side,
732 auto row_type, auto col_type) {
733 return add_bmc.get<0>().find(boost::make_tuple(
734 row_name, col_name, row_type, col_type, row_side, col_side));
735 };
736
737 auto set_block = [&](BlockMatContainer &add_bmc, auto &row_name,
738 auto &col_name, auto row_side, auto col_side,
739 auto row_type, auto col_type, const auto &m,
740 const auto &row_ind, const auto &col_ind) {
741 auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
742 row_type, col_type);
743 if (it != add_bmc.get<0>().end()) {
744 it->setMat(m);
745 it->setInd(row_ind, col_ind);
746 it->setSetAtElement();
747 return it;
748 } else {
749 auto p_it = add_bmc.insert(BlockMatData(row_name, col_name, row_type,
750 col_type, row_side, col_side,
751 m, row_ind, col_ind));
752 return p_it.first;
753 }
754 };
755
756 auto add_block = [&](BlockMatContainer &add_bmc, auto &row_name,
757 auto &col_name, auto row_side, auto col_side,
758 auto row_type, auto col_type, const auto &m,
759 const auto &row_ind, const auto &col_ind) {
760 auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
761 row_type, col_type);
762 if (it != add_bmc.get<0>().end()) {
763 it->addMat(m);
764 it->setInd(row_ind, col_ind);
765 it->setSetAtElement();
766 return it;
767 } else {
768 auto p_it = add_bmc.insert(BlockMatData(row_name, col_name, row_type,
769 col_type, row_side, col_side,
770 m, row_ind, col_ind));
771 return p_it.first;
772 }
773 };
774
775 auto assemble_block = [&](auto &bit, Mat S) {
777 const VectorInt &rind = bit.rowInd;
778 const VectorInt &cind = bit.colInd;
779 const MatrixDouble &m = bit.M;
780
781 CHKERR MatSetValues(S, rind.size(), &*rind.begin(), cind.size(),
782 &*cind.begin(), &*m.data().begin(), ADD_VALUES);
783
785 };
786
787 auto create_block_schur_bc = [&](BlockMatContainer &bmc,
788 BlockMatContainer &add_bmc,
789 std::string field, AO ao) {
791 for (auto &bit : add_bmc) {
792 bit.unSetAtElement();
793 bit.clearMat();
794 }
795
796 for (auto &bit : bmc) {
797 if (bit.setAtElement && bit.rowField != field &&
798 bit.colField != field) {
799 VectorInt rind = bit.rowInd;
800 VectorInt cind = bit.colInd;
801 const MatrixDouble &m = bit.M;
802
803 // for (int i = 0; i != m.size2() * m.size1(); i++)
804 // if (abs(m.data()[i]) > 1e-8) {
805 // string wait;
806 // cout << "what" << endl;
807 // }
808
809 if (ao) {
810 CHKERR AOApplicationToPetsc(ao, rind.size(), &*rind.begin());
811 CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
812 }
813 auto it =
814 set_block(add_bmc, bit.rowField, bit.colField, bit.rowSide,
815 bit.colSide, bit.rowType, bit.colType, m, rind, cind);
816 }
817 }
818
819 // for (auto &bit_col : bmc) {
820 // if (bit_col.setAtElement && bit_col.rowField == field &&
821 // bit_col.colField != field) {
822 // const MatrixDouble &cm = bit_col.M;
823 // VectorInt cind = bit_col.colInd;
824 // // invMat.resize(inv.size1(), cm.size2(), false);
825 // // noalias(invMat) = prod(inv, cm);
826 // noalias(invMat) = cm;
827 // if (ao)
828 // CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
829 // for (auto &bit_row : bmc) {
830 // if (bit_row.setAtElement && bit_row.rowField != field &&
831 // bit_row.colField == field) {
832 // const MatrixDouble &rm = bit_row.M;
833 // VectorInt rind = bit_row.rowInd;
834 // K.resize(rind.size(), cind.size(), false);
835 // // noalias(K) = prod(rm, invMat);
836 // noalias(K) = rm;
837 // if (ao)
838 // CHKERR AOApplicationToPetsc(ao, rind.size(), &*rind.begin());
839 // auto it =
840 // add_block(add_bmc, bit_row.rowField, bit_col.colField,
841 // bit_row.rowSide, bit_col.colSide, bit_row.rowType,
842 // bit_col.colType, K, rind, cind);
843 // }
844 // }
845 // }
846 // }
847
849 };
850
851 auto assemble_schur_bc = [&](BlockMatContainer &add_bmc, Mat S,
852 bool debug = false) {
854 for (auto &bit : add_bmc) {
855 if (bit.setAtElement)
856 CHKERR assemble_block(bit, S);
857 }
858 if (debug) {
859 for (auto &bit : add_bmc) {
860 if (bit.setAtElement) {
861 std::cerr << "assemble: " << bit.rowField << " " << bit.colField
862 << endl;
863 // std::cerr << bit.M << endl;
864 }
865 }
866 std::cerr << std::endl;
867 }
869 };
870
871 auto precondition_schur = [&](BlockMatContainer &bmc,
872 BlockMatContainer &add_bmc,
873 const std::string field,
874 const MatrixDouble &diag_mat,
875 const double eps) {
877
878 for (auto &bit : add_bmc) {
879 bit.unSetAtElement();
880 bit.clearMat();
881 }
882
883 for (auto &bit : bmc) {
884 if (bit.setAtElement) {
885 if (bit.rowField != field || bit.colField != field)
886 auto it = set_block(add_bmc, bit.rowField, bit.colField,
887 bit.rowSide, bit.colSide, bit.rowType,
888 bit.colType, bit.M, bit.rowInd, bit.colInd);
889 }
890 }
891
892 auto bit =
893 bmc.get<1>().find(boost::make_tuple(field, field, row_type, row_type));
894 if (bit->setAtElement && bit != bmc.get<1>().end()) {
895 auto it = set_block(add_bmc, bit->rowField, bit->colField,
896 bit->rowSide, bit->colSide, bit->rowType,
897 bit->colType, bit->M, bit->rowInd, bit->colInd);
898 MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
899 m += eps * diag_mat;
900 } else {
901 auto row_it = bmc.get<3>().lower_bound(field);
902 for (; row_it != bmc.get<3>().upper_bound(field); ++row_it) {
903 if (row_it->setAtElement) {
904 auto it = set_block(add_bmc, field, field, 0, 0, row_type, row_type,
905 diag_mat, row_it->rowInd, row_it->rowInd);
906 MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
907 m *= eps;
908 break;
909 }
910 }
911 if (row_it == bmc.get<3>().end())
912 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
913 "row field not found %s", field.c_str());
914 }
915
917 };
918
919 const bool debug = false;
920 auto &mat_cont = *commonDataPtr->blockMatContainerPtr;
921 auto &mass_mat_tensor = commonDataPtr->massMatTensor;
922
923 if (SEpEp) {
924 CHKERR create_block_schur_bc(mat_cont, blockMat["BC"], "EP", aoSEpEp);
925 // if (getPtrFE()->mField.check_field("SIGMA") && false) {
926 // constexpr double eps = 1e-8;
927 // CHKERR precondition_schur(blockMat["BC"], blockMat["precBC"],
928 // "SIGMA",
929 // mass_mat_tensor, eps);
930 // CHKERR assemble_schur_bc(blockMat["precBC"], SEpEp, debug);
931 // blockMat["BC"] = blockMat["precBC"];
932 // // CHKERR assemble_schur_bc(blockMat["BC"], SEpEp, debug);
933 // } else
934 CHKERR assemble_schur_bc(blockMat["BC"], SEpEp, debug);
935
936 if (STauTau) {
937 CHKERR create_block_schur_bc(blockMat["BC"], blockMat["TAUTAU"],
938 "TAU", aoSTauTau);
939 CHKERR assemble_schur_bc(blockMat["TAUTAU"], STauTau, debug);
940 }
941 }
942
943 blockMat.clear();
944 }
945 }
946
948}
949
951 EntData &data) {
953 if (row_type != MBTET && row_type != MBHEX)
955 if (!commonDataPtr->blockMatContainerPtr)
957 auto &bmc = *commonDataPtr->blockMatContainerPtr;
958 for (auto &bit : bmc)
959 bit.unSetAtElement();
960 // bmc.clear();
962}
963
965 EntData &data) {
967 if (row_type != MBTET && row_type != MBHEX)
969 if (commonDataPtr->blockMatContainerPtr) {
970
971 auto SEpEp = commonDataPtr->SEpEp;
972 auto aoSEpEp = commonDataPtr->aoSEpEp;
973 auto STauTau = commonDataPtr->STauTau;
974 auto aoSTauTau = commonDataPtr->aoSTauTau;
975
976 if (SEpEp) {
977
978 auto find_block = [&](BlockMatContainer &add_bmc, auto &row_name,
979 auto &col_name, auto row_side, auto col_side,
980 auto row_type, auto col_type) {
981 return add_bmc.get<0>().find(boost::make_tuple(
982 row_name, col_name, row_type, col_type, row_side, col_side));
983 };
984
985 auto set_block = [&](BlockMatContainer &add_bmc, auto &row_name,
986 auto &col_name, auto row_side, auto col_side,
987 auto row_type, auto col_type, const auto &m,
988 const auto &row_ind, const auto &col_ind) {
989 auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
990 row_type, col_type);
991 if (it != add_bmc.get<0>().end()) {
992 it->setMat(m);
993 it->setInd(row_ind, col_ind);
994 it->setSetAtElement();
995 return it;
996 } else {
997 auto p_it = add_bmc.insert(BlockMatData(row_name, col_name, row_type,
998 col_type, row_side, col_side,
999 m, row_ind, col_ind));
1000 return p_it.first;
1001 }
1002 };
1003
1004 auto add_block = [&](BlockMatContainer &add_bmc, auto &row_name,
1005 auto &col_name, auto row_side, auto col_side,
1006 auto row_type, auto col_type, const auto &m,
1007 const auto &row_ind, const auto &col_ind) {
1008 auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
1009 row_type, col_type);
1010 if (it != add_bmc.get<0>().end()) {
1011 it->addMat(m);
1012 it->setInd(row_ind, col_ind);
1013 it->setSetAtElement();
1014 return it;
1015 } else {
1016 auto p_it = add_bmc.insert(BlockMatData(row_name, col_name, row_type,
1017 col_type, row_side, col_side,
1018 m, row_ind, col_ind));
1019 return p_it.first;
1020 }
1021 };
1022
1023 auto assemble_block = [&](auto &bit, Mat S) {
1025 const VectorInt &rind = bit.rowInd;
1026 const VectorInt &cind = bit.colInd;
1027 const MatrixDouble &m = bit.M;
1028
1029 CHKERR MatSetValues(S, rind.size(), &*rind.begin(), cind.size(),
1030 &*cind.begin(), &*m.data().begin(), ADD_VALUES);
1031
1033 };
1034
1035 auto invert_symm_mat = [&](MatrixDouble &m, auto &inv) {
1037 const int nb = m.size1();
1038
1039 inv.resize(nb, nb, false);
1040 inv.clear();
1041 for (int cc = 0; cc != nb; ++cc)
1042 inv(cc, cc) = -1;
1043
1044 iPIV.resize(nb, false);
1045 lapackWork.resize(nb * nb, false);
1046 const auto info = lapack_dsysv('L', nb, nb, &*m.data().begin(), nb,
1047 &*iPIV.begin(), &*inv.data().begin(), nb,
1048 &*lapackWork.begin(), nb * nb);
1049 // if (info != 0)
1050 // SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1051 // "Can not invert matrix info = %d", info);
1052
1054 };
1055
1056 auto invert_symm_schur = [&](BlockMatContainer &bmc, std::string field,
1057 auto &inv) {
1059
1060 auto bit =
1061 bmc.get<1>().find(boost::make_tuple(field, field, row_type, row_type));
1062 if (bit != bmc.get<1>().end()) {
1063
1064 auto &m = *const_cast<MatrixDouble *>(&(bit->M));
1065 CHKERR invert_symm_mat(m, inv);
1066
1067 } else
1068 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1069 "%s matrix not found", field.c_str());
1070
1072 };
1073
1074 // auto invert_nonsymm_mat = [&](MatrixDouble &m, auto &inv) {
1075 // MoFEMFunctionBeginHot;
1076
1077 // const int nb = m.size1();
1078
1079 // MatrixDouble trans_m = trans(m);
1080 // MatrixDouble trans_inv;
1081 // trans_inv.resize(nb, nb, false);
1082 // trans_inv.clear();
1083 // for (int c = 0; c != nb; ++c)
1084 // trans_inv(c, c) = -1;
1085
1086 // iPIV.resize(nb, false);
1087 // const auto info =
1088 // lapack_dgesv(nb, nb, &*trans_m.data().begin(), nb,
1089 // &*iPIV.begin(),
1090 // &*trans_inv.data().begin(), nb);
1091 // if (info != 0)
1092 // SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1093 // "Can not invert matrix info = %d", info);
1094
1095 // inv.resize(nb, nb, false);
1096 // noalias(inv) = trans(trans_inv);
1097
1098 // MoFEMFunctionReturnHot(0);
1099 // };
1100
1101 // auto invert_nonsymm_schur = [&](BlockMatContainer &bmc,
1102 // std::string field, auto &inv,
1103 // const bool debug = false) {
1104 // MoFEMFunctionBeginHot;
1105
1106 // auto bit =
1107 // bmc.get<1>().find(boost::make_tuple(field, field, row_type, row_type));
1108 // if (bit != bmc.get<1>().end()) {
1109
1110 // auto &m = *const_cast<MatrixDouble *>(&(bit->M));
1111 // CHKERR invert_nonsymm_mat(m, inv);
1112
1113 // if (debug) {
1114 // std::cerr << prod(m, inv) << endl;
1115 // std::cerr << endl;
1116 // }
1117
1118 // } else
1119 // SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1120 // "%s matrix not found", field.c_str());
1121
1122 // MoFEMFunctionReturnHot(0);
1123 // };
1124
1125 auto create_block_schur = [&](BlockMatContainer &bmc,
1126 BlockMatContainer &add_bmc,
1127 std::string field, AO ao,
1128 MatrixDouble &inv) {
1130 // AOView(ao, PETSC_VIEWER_STDOUT_SELF);
1131 for (auto &bit : add_bmc) {
1132 bit.unSetAtElement();
1133 bit.clearMat();
1134 }
1135
1136 for (auto &bit : bmc) {
1137 if (bit.setAtElement && bit.rowField != field &&
1138 bit.colField != field) {
1139 VectorInt rind = bit.rowInd;
1140 VectorInt cind = bit.colInd;
1141 const MatrixDouble &m = bit.M;
1142 if (ao) {
1143 CHKERR AOApplicationToPetsc(ao, rind.size(), &*rind.begin());
1144 CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
1145 }
1146 auto it =
1147 set_block(add_bmc, bit.rowField, bit.colField, bit.rowSide,
1148 bit.colSide, bit.rowType, bit.colType, m, rind, cind);
1149 }
1150 }
1151
1152 for (auto &bit_col : bmc) {
1153 if (bit_col.setAtElement && bit_col.rowField == field &&
1154 bit_col.colField != field) {
1155 const MatrixDouble &cm = bit_col.M;
1156 VectorInt cind = bit_col.colInd;
1157 invMat.resize(inv.size1(), cm.size2(), false);
1158 noalias(invMat) = prod(inv, cm);
1159 if (ao)
1160 CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
1161 for (auto &bit_row : bmc) {
1162 if (bit_row.setAtElement && bit_row.rowField != field &&
1163 bit_row.colField == field) {
1164 const MatrixDouble &rm = bit_row.M;
1165 VectorInt rind = bit_row.rowInd;
1166 K.resize(rind.size(), cind.size(), false);
1167 noalias(K) = prod(rm, invMat);
1168 if (ao)
1169 CHKERR AOApplicationToPetsc(ao, rind.size(), &*rind.begin());
1170 auto it =
1171 add_block(add_bmc, bit_row.rowField, bit_col.colField,
1172 bit_row.rowSide, bit_col.colSide, bit_row.rowType,
1173 bit_col.colType, K, rind, cind);
1174 }
1175 }
1176 }
1177 }
1178
1180 };
1181
1182 auto assemble_schur = [&](BlockMatContainer &add_bmc, Mat S,
1183 bool debug = false) {
1185 if (debug) {
1186 for (auto &bit : add_bmc) {
1187 if (bit.setAtElement) {
1188 std::cerr << "assemble vol: " << bit.rowField << " "
1189 << bit.colField << endl;
1190 // std::cerr << bit.M << endl;
1191 }
1192 }
1193 std::cerr << std::endl;
1194 }
1195 for (auto &bit : add_bmc) {
1196 if (bit.setAtElement)
1197 CHKERR assemble_block(bit, S);
1198 }
1200 };
1201
1202 auto precondition_schur = [&](BlockMatContainer &bmc,
1203 BlockMatContainer &add_bmc,
1204 const std::string field,
1205 const MatrixDouble &diag_mat,
1206 const double eps) {
1208
1209 for (auto &bit : add_bmc) {
1210 bit.unSetAtElement();
1211 bit.clearMat();
1212 }
1213
1214 for (auto &bit : bmc) {
1215 if (bit.setAtElement) {
1216 if (bit.rowField != field || bit.colField != field)
1217 auto it = set_block(add_bmc, bit.rowField, bit.colField,
1218 bit.rowSide, bit.colSide, bit.rowType,
1219 bit.colType, bit.M, bit.rowInd, bit.colInd);
1220 }
1221 }
1222
1223 auto bit =
1224 bmc.get<1>().find(boost::make_tuple(field, field, row_type, row_type));
1225 if (bit->setAtElement && bit != bmc.get<1>().end()) {
1226 auto it = set_block(add_bmc, bit->rowField, bit->colField,
1227 bit->rowSide, bit->colSide, bit->rowType,
1228 bit->colType, bit->M, bit->rowInd, bit->colInd);
1229 MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
1230 m += eps * diag_mat;
1231 } else {
1232 auto row_it = bmc.get<3>().lower_bound(field);
1233 for (; row_it != bmc.get<3>().upper_bound(field); ++row_it) {
1234 if (row_it->setAtElement) {
1235 auto it = set_block(add_bmc, field, field, 0, 0, row_type, row_type,
1236 diag_mat, row_it->rowInd, row_it->rowInd);
1237 MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
1238 m *= eps;
1239 break;
1240 }
1241 }
1242 if (row_it == bmc.get<3>().end())
1243 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1244 "row field not found %s", field.c_str());
1245 }
1246
1248 };
1249
1250 const bool debug = false;
1251 auto &mat_cont = *commonDataPtr->blockMatContainerPtr;
1252
1253 CHKERR invert_symm_schur(mat_cont, "EP",
1254 invBlockMat["EPEP"]);
1255 CHKERR create_block_schur(mat_cont,
1256 blockMat["EPEP"], "EP", aoSEpEp,
1257 invBlockMat["EPEP"]);
1258
1259 auto &mass_mat_scalar = commonDataPtr->massMatScalar;
1260 constexpr double eps = 1e-8; //FIXME: TODO: figure out this later
1261 CHKERR precondition_schur(blockMat["EPEP"], blockMat["precEPEP"], "TAU",
1262 mass_mat_scalar, eps);
1263 // if (getPtrFE()->mField.check_field("SIGMA") && false) {
1264 // auto &mass_mat_tensor = commonDataPtr->massMatTensor;
1265 // CHKERR precondition_schur(blockMat["precEPEP"], blockMat["precSIGMA"],
1266 // "SIGMA", mass_mat_tensor, eps);
1267 // CHKERR assemble_schur(blockMat["precSIGMA"], SEpEp, debug);
1268 // blockMat["precEPEP"] = blockMat["precSIGMA"];
1269 // } else
1270 CHKERR assemble_schur(blockMat["precEPEP"], SEpEp, false);
1271
1272 if(STauTau){
1273 CHKERR invert_symm_schur(blockMat["precEPEP"], "TAU",
1274 invBlockMat["TAUTAU"]);
1275 CHKERR create_block_schur(blockMat["precEPEP"], blockMat["TAUTAU"], "TAU",
1276 aoSTauTau, invBlockMat["TAUTAU"]);
1277 CHKERR assemble_schur(blockMat["TAUTAU"], STauTau);
1278 }
1279
1280 blockMat.clear();
1281
1282 }
1283 }
1284
1286}
1287
1288template struct OpInternalForceRhs<true, true>;
1289template struct OpInternalForceRhs<false, false>;
1290template struct OpInternalForceRhs<true, false>;
1291template struct OpInternalForceRhs<false, true>;
1292
1293template struct OpPostProcElastic<true>;
1294template struct OpPostProcElastic<false>;
1295
1296template struct OpCalculateLogStressTangent<true>;
1298
1299} // namespace OpElasticTools
static Index< 'o', 3 > o
static Index< 'J', 3 > J
static Index< 'p', 3 > p
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_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#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
static const bool debug
#define MU(E, NU)
Definition: fem_tools.h:23
#define LAMBDA(E, NU)
Definition: fem_tools.h:22
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
constexpr double lambda
constexpr auto t_kd
MoFEMErrorCode VecSetValues< EssentialBcStorage >(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Set values to vector in operator.
auto bit
set bit
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
static __CLPK_integer lapack_dsysv(char uplo, __CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition: lapack_wrap.h:220
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
std::map< int, BlockParamData > mat_blocks
BlockParamData * cache
Tensors class implemented by Walter Landry.
Definition: FTensor.hpp:51
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasVector< int > VectorInt
Definition: Types.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:1323
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
multi_index_container< BlockMatData, indexed_by< ordered_unique< composite_key< BlockMatData, member< BlockMatData, std::string, &BlockMatData::rowField >, member< BlockMatData, std::string, &BlockMatData::colField >, member< BlockMatData, EntityType, &BlockMatData::rowType >, member< BlockMatData, EntityType, &BlockMatData::colType >, member< BlockMatData, int, &BlockMatData::rowSide >, member< BlockMatData, int, &BlockMatData::colSide > > >, ordered_non_unique< composite_key< BlockMatData, member< BlockMatData, std::string, &BlockMatData::rowField >, member< BlockMatData, std::string, &BlockMatData::colField >, member< BlockMatData, EntityType, &BlockMatData::rowType >, member< BlockMatData, EntityType, &BlockMatData::colType > > >, ordered_non_unique< composite_key< BlockMatData, member< BlockMatData, std::string, &BlockMatData::rowField >, member< BlockMatData, std::string, &BlockMatData::colField > > >, ordered_non_unique< member< BlockMatData, std::string, &BlockMatData::rowField > >, ordered_non_unique< member< BlockMatData, std::string, &BlockMatData::colField > > > > BlockMatContainer
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1312
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_lapack(const Tensor2_symmetric< T1, 3 > &X, Tensor1< T2, 3 > &lambda, Tensor2< T3, 3, 3 > &eig_vec)
Tensor4< double, 3, 3, 3, 3 > getDTensorFunction(CommonData &common_data, Tensor2< T, 3, 3 > &F, Tensor2< double, 3, 3 > &t_stress)
auto to_non_symm(T &symm)
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)
constexpr double t
plate stiffness
Definition: plate.cpp:59
constexpr auto field_name
constexpr double mu
chrono::high_resolution_clock::time_point stop
chrono::high_resolution_clock::time_point start
virtual moab::Interface & get_moab()=0
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
const VectorInt & getIndices() const
Get global indices of dofs on entity.
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateLogStressTangent(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > mat_D_ptr)
OpCalculateNeoHookeStressTangent(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode iNtegrate(EntData &data)
[Internal force]
MoFEMErrorCode aSsemble(EntData &data)
OpInternalForceRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
[Calculate stress]
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Calculate strain]
boost::shared_ptr< CommonData > commonDataPtr
OpLogStrain(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
OpPostProcElastic(const std::string field_name, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, boost::shared_ptr< CommonData > common_data_ptr)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Postprocessing]
OpSaveReactionForces(MoFEM::Interface &m_field, const std::string field_name, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, boost::shared_ptr< CommonData > common_data_ptr)
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
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
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
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
map< std::string, MatrixDouble > invBlockMat
map< std::string, BlockMatContainer > blockMat
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.
OpSetMaterialBlockBoundary(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr, std::map< EntityHandle, int > &map_block_tets)
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
OpSetMaterialBlock(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr, std::map< EntityHandle, int > &map_block_tets)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Calculate strain]
OpStrain(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
OpStress(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > mat_D_Ptr)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Calculate stress]
boost::shared_ptr< MatrixDouble > matDPtr