35 chrono::high_resolution_clock::time_point
start;
36 chrono::high_resolution_clock::time_point
stop;
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;
48 boost::shared_ptr<CommonData> common_data_ptr)
50 commonDataPtr(common_data_ptr) {
60 auto t_grad = getFTensor2FromMat<3, 3>(*(
commonDataPtr->mGradPtr));
61 auto t_strain = getFTensor2SymmetricFromMat<3>(*(
commonDataPtr->mStrainPtr));
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;
77 boost::shared_ptr<CommonData> common_data_ptr)
79 commonDataPtr(common_data_ptr) {
89 commonDataPtr->mDeformationGradient->resize(9, nb_gauss_pts,
false);
92 commonDataPtr->matEigenVector->resize(9, nb_gauss_pts,
false);
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));
112 auto t_strain = getFTensor2SymmetricFromMat<3>(*(
commonDataPtr->mStrainPtr));
114 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
137 t_strain(
i,
j) = 0.5 * lnC(
i,
j);
151 boost::shared_ptr<CommonData> common_data_ptr, boost::shared_ptr<MatrixDouble> mat_D_Ptr)
153 commonDataPtr(common_data_ptr), matDPtr(mat_D_Ptr) {
161 const size_t nb_gauss_pts =
commonDataPtr->mStrainPtr->size2();
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));
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);
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) {
185template <
bool TANGENT>
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));
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);
206 dE_dF.resize(81, nb_gauss_pts,
false);
207 auto t_dE_dF = getFTensor4FromMat<3, 3, 3, 3>(dE_dF);
214 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
217 auto dlogC_dC = calculate_lagrangian_moduli_TL<TANGENT>(t_stress, t_eigen,
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;
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);
228 P_D_P_plus_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;
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));
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);
250 if (
false && TANGENT) {
251 if (getTStime() > 0.2) {
254 cout <<
"NEW IMPLEMENTATION" << endl;
256 for (
int ii = 0; ii != 1000; ii++)
258 C, t_stress, t_eigen, t_eigen_vec, TL3_lukasz);
262 cout <<
"NEW IMPLEMENTATION" << endl;
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);
269 cout <<
"(OLD) FINITE DIFFERENCE" << endl;
272 for (
int ii = 0; ii != 1000; ii++)
276 cout <<
"EIGENVALUES (LAPACK)" << endl;
278 for (
int ii = 0; ii != 1000; ii++)
302 if (getTStime() > 0.2)
320 const std::string
field_name, boost::shared_ptr<CommonData> common_data_ptr)
322 commonDataPtr(common_data_ptr) {
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));
339 dP_dF.resize(81, nb_gauss_pts,
false);
340 auto t_dP_dF = getFTensor4FromMat<3, 3, 3, 3>(dP_dF);
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;
348 const double lambda =
LAMBDA((*cache).young_modulus, (*cache).poisson);
349 const double mu =
MU((*cache).young_modulus, (*cache).poisson);
351 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
355 const double log_det = log(t_detF);
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);
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;
376 const double var =
lambda * log_det -
mu;
377 t_piola(
i,
J) =
mu *
F(
i,
J) + var * invF(
J,
i);
391template <
bool LOGSTRAIN,
bool REACTIONS>
393 const std::string
field_name, boost::shared_ptr<CommonData> common_data_ptr)
395 commonDataPtr(common_data_ptr) {}
398template <
bool LOGSTRAIN,
bool REACTIONS>
403 const size_t nb_dofs = data.
getIndices().size();
406 const size_t nb_base_functions = data.
getN().size2();
407 if (3 * nb_base_functions < nb_dofs)
410 "Number of DOFs is larger than number of base functions on entity");
412 const size_t nb_gauss_pts = data.
getN().size1();
414 auto t_w = getFTensor0IntegrationWeight();
416 getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStressPtr));
417 auto t_piola = getFTensor2FromMat<3, 3>(*(commonDataPtr->mPiolaStressPtr));
420 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
422 double alpha = getMeasure() * t_w;
423 auto t_nf = getFTensor1FromArray<3, 3>(locF);
427 out_stress(
i,
j) = t_piola(
i,
j);
432 for (; bb != nb_dofs / 3; ++bb) {
433 t_nf(
i) += alpha * t_diff_base(
j) * out_stress(
i,
j);
437 for (; bb < nb_base_functions; ++bb)
447template <
bool LOGSTRAIN,
bool REACTIONS>
452 auto react_ents = commonDataPtr->reactionEnts;
455 auto &vec = commonDataPtr->reactionsVec;
457 CHKERR VecGetArray(vec, &react_ptr);
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,
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;
475 CHKERR m_field.get_moab().tag_set_data(commonDataPtr->reactionTag, &ent,
478 CHKERR VecRestoreArray(vec, &react_ptr);
481 this->getKSPf(), data, &*this->locF.data().begin(), ADD_VALUES);
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) {
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) {
524 std::vector<EntityHandle> adj_ents;
527 CHKERR moab.get_adjacencies(&fe_ent, 1, 3,
false, adj_ents);
538 moab::Interface &post_proc_mesh, std::vector<EntityHandle> &map_gauss_pts,
539 boost::shared_ptr<CommonData> common_data_ptr)
541 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
542 commonDataPtr(common_data_ptr) {
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);
556 MB_TAG_CREAT | MB_TAG_SPARSE,
561 auto th_react = get_tag(
"REACTIONS", 3);
562 size_t nb_gauss_pts = data.
getN().size1();
565 const size_t rank = data.
getFieldDofs()[0]->getNbOfCoeffs();
569 auto ent = dof->getEnt();
571 const size_t idx = dof->getDofCoeffIdx();
573 double tag_val[rank];
576 vAlues(
dd++) = tag_val[idx];
582 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
583 for (
int rr = 0; rr != rank; rr++) {
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;
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) {
606template <
bool LOGSTRAIN>
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);
615 CHKERR postProcMesh.tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE,
th,
616 MB_TAG_CREAT | MB_TAG_SPARSE,
625 for (
size_t r = 0; r != 3; ++r)
626 for (
size_t c = 0;
c != 3; ++
c)
638 return postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1,
639 &*mat.data().begin());
642 string strain_tag_name =
"STRAIN";
643 string stress_tag_name =
"CAUCHY_STRESS";
645 strain_tag_name =
"LOG_STRAIN";
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);
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));
661 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
663 t_stress(
i,
j) /= (*cache).young_modulus_inv;
667 t_piola(
i,
j) /= (*cache).young_modulus_inv;
672 t_stress(
i,
j) = (1. / detF) * (t_piola(
i,
k) ^
F(
j,
k));
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))));
683 CHKERR set_tag(th_mises, gg, set_scalar(mises));
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));
701 if (row_type != MBTRI && row_type != MBQUAD)
708 for (
auto &
bit : bmc)
709 bit.unSetAtElement();
716 if (row_type != MBTRI && row_type != MBQUAD)
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));
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,
743 if (it != add_bmc.get<0>().end()) {
745 it->setInd(row_ind, col_ind);
746 it->setSetAtElement();
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));
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,
762 if (it != add_bmc.get<0>().end()) {
764 it->setInd(row_ind, col_ind);
765 it->setSetAtElement();
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));
775 auto assemble_block = [&](
auto &
bit, Mat S) {
782 &*cind.begin(), &*
m.data().begin(), ADD_VALUES);
789 std::string field, AO ao) {
791 for (
auto &
bit : add_bmc) {
792 bit.unSetAtElement();
796 for (
auto &
bit : bmc) {
797 if (
bit.setAtElement &&
bit.rowField != field &&
798 bit.colField != field) {
810 CHKERR AOApplicationToPetsc(ao, rind.size(), &*rind.begin());
811 CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
814 set_block(add_bmc,
bit.rowField,
bit.colField,
bit.rowSide,
815 bit.colSide,
bit.rowType,
bit.colType,
m, rind, cind);
852 bool debug =
false) {
854 for (
auto &
bit : add_bmc) {
855 if (
bit.setAtElement)
859 for (
auto &
bit : add_bmc) {
860 if (
bit.setAtElement) {
861 std::cerr <<
"assemble: " <<
bit.rowField <<
" " <<
bit.colField
866 std::cerr << std::endl;
873 const std::string field,
878 for (
auto &
bit : add_bmc) {
879 bit.unSetAtElement();
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,
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,
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);
911 if (row_it == bmc.get<3>().end())
913 "row field not found %s", field.c_str());
919 const bool debug =
false;
924 CHKERR create_block_schur_bc(mat_cont,
blockMat[
"BC"],
"EP", aoSEpEp);
953 if (row_type != MBTET && row_type != MBHEX)
958 for (
auto &
bit : bmc)
959 bit.unSetAtElement();
967 if (row_type != MBTET && row_type != MBHEX)
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));
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,
991 if (it != add_bmc.get<0>().end()) {
993 it->setInd(row_ind, col_ind);
994 it->setSetAtElement();
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));
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()) {
1012 it->setInd(row_ind, col_ind);
1013 it->setSetAtElement();
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));
1023 auto assemble_block = [&](
auto &
bit, Mat S) {
1030 &*cind.begin(), &*
m.data().begin(), ADD_VALUES);
1037 const int nb =
m.size1();
1039 inv.resize(nb, nb,
false);
1041 for (
int cc = 0; cc != nb; ++cc)
1044 iPIV.resize(nb,
false);
1046 const auto info =
lapack_dsysv(
'L', nb, nb, &*
m.data().begin(), nb,
1047 &*
iPIV.begin(), &*inv.data().begin(), nb,
1061 bmc.get<1>().find(boost::make_tuple(field, field, row_type, row_type));
1062 if (
bit != bmc.get<1>().end()) {
1065 CHKERR invert_symm_mat(
m, inv);
1069 "%s matrix not found", field.c_str());
1127 std::string field, AO ao,
1131 for (
auto &
bit : add_bmc) {
1132 bit.unSetAtElement();
1136 for (
auto &
bit : bmc) {
1137 if (
bit.setAtElement &&
bit.rowField != field &&
1138 bit.colField != field) {
1143 CHKERR AOApplicationToPetsc(ao, rind.size(), &*rind.begin());
1144 CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
1147 set_block(add_bmc,
bit.rowField,
bit.colField,
bit.rowSide,
1148 bit.colSide,
bit.rowType,
bit.colType,
m, rind, cind);
1152 for (
auto &bit_col : bmc) {
1153 if (bit_col.setAtElement && bit_col.rowField == field &&
1154 bit_col.colField != field) {
1157 invMat.resize(inv.size1(), cm.size2(),
false);
1158 noalias(
invMat) = prod(inv, cm);
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) {
1166 K.resize(rind.size(), cind.size(),
false);
1167 noalias(
K) = prod(rm,
invMat);
1169 CHKERR AOApplicationToPetsc(ao, rind.size(), &*rind.begin());
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);
1183 bool debug =
false) {
1186 for (
auto &
bit : add_bmc) {
1187 if (
bit.setAtElement) {
1188 std::cerr <<
"assemble vol: " <<
bit.rowField <<
" "
1189 <<
bit.colField << endl;
1193 std::cerr << std::endl;
1195 for (
auto &
bit : add_bmc) {
1196 if (
bit.setAtElement)
1204 const std::string field,
1209 for (
auto &
bit : add_bmc) {
1210 bit.unSetAtElement();
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,
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,
1230 m +=
eps * diag_mat;
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);
1242 if (row_it == bmc.get<3>().end())
1244 "row field not found %s", field.c_str());
1250 const bool debug =
false;
1253 CHKERR invert_symm_schur(mat_cont,
"EP",
1255 CHKERR create_block_schur(mat_cont,
1260 constexpr double eps = 1e-8;
1262 mass_mat_scalar,
eps);
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
constexpr double lambda
surface tension
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)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
std::map< int, BlockParamData > mat_blocks
Tensors class implemented by Walter Landry.
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)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasVector< int > VectorInt
implementation of Data Operators for Forces and Sources
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.
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.
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.
constexpr double t
plate stiffness
constexpr auto field_name
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.
default operator for TRI element
const TSMethod::TSContext getTSCtx() const
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
ForcesAndSourcesCore * getPtrFE() const
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]