18#ifndef __HOOKE_ELEMENT_HPP
19#define __HOOKE_ELEMENT_HPP
21#ifndef __BASICFINITEELEMENTS_HPP__
25#ifndef __NONLINEAR_ELASTIC_HPP
56#ifndef __CONVECTIVE_MASS_ELEMENT_HPP
77 VolumeElementForcesAndSourcesCore::UserDataOperator;
82 boost::shared_ptr<MatrixDouble>
hMat;
83 boost::shared_ptr<MatrixDouble>
FMat;
85 boost::shared_ptr<MatrixDouble>
HMat;
98 smallStrainMat = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
99 hMat = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
100 FMat = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
102 HMat = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
103 detHVec = boost::shared_ptr<VectorDouble>(
new VectorDouble());
104 invHMat = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
106 cauchyStressMat = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
107 stiffnessMat = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
108 energyVec = boost::shared_ptr<VectorDouble>(
new VectorDouble());
109 eshelbyStressMat = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
111 eshelbyStress_dx = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
118 template <
bool D = false>
122 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
133 const std::string col_field,
134 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
142#define MAT_TO_DDG(SM) \
143 &(*SM)(0, 0), &(*SM)(1, 0), &(*SM)(2, 0), &(*SM)(3, 0), &(*SM)(4, 0), \
144 &(*SM)(5, 0), &(*SM)(6, 0), &(*SM)(7, 0), &(*SM)(8, 0), &(*SM)(9, 0), \
145 &(*SM)(10, 0), &(*SM)(11, 0), &(*SM)(12, 0), &(*SM)(13, 0), \
146 &(*SM)(14, 0), &(*SM)(15, 0), &(*SM)(16, 0), &(*SM)(17, 0), \
147 &(*SM)(18, 0), &(*SM)(19, 0), &(*SM)(20, 0), &(*SM)(21, 0), \
148 &(*SM)(22, 0), &(*SM)(23, 0), &(*SM)(24, 0), &(*SM)(25, 0), \
149 &(*SM)(26, 0), &(*SM)(27, 0), &(*SM)(28, 0), &(*SM)(29, 0), \
150 &(*SM)(30, 0), &(*SM)(31, 0), &(*SM)(32, 0), &(*SM)(33, 0), \
151 &(*SM)(34, 0), &(*SM)(35, 0)
156 boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
167 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
180 const std::string row_field,
const std::string col_field,
181 boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
193 const std::string row_field,
const std::string col_field,
194 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
195 boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
200 boost::shared_ptr<map<int, BlockData>>
223 const std::string col_field,
BlockData &data,
225 boost::shared_ptr<DataAtIntegrationPts> &common_data,
228 row_field, col_field, OPROWCOL, symm),
229 commonData(common_data), dAta(data), massData(mass_data) {}
237 auto get_tensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
239 &
m(3 * r + 0, 3 *
c + 0), &
m(3 * r + 0, 3 *
c + 1),
240 &
m(3 * r + 0, 3 *
c + 2), &
m(3 * r + 1, 3 *
c + 0),
241 &
m(3 * r + 1, 3 *
c + 1), &
m(3 * r + 1, 3 *
c + 2),
242 &
m(3 * r + 2, 3 *
c + 0), &
m(3 * r + 2, 3 *
c + 1),
243 &
m(3 * r + 2, 3 *
c + 2));
246 const int row_nb_dofs = row_data.
getIndices().size();
249 const int col_nb_dofs = col_data.
getIndices().size();
252 if (dAta.tEts.find(getFEEntityHandle()) == dAta.tEts.end()) {
255 if (massData.
tEts.find(getFEEntityHandle()) == massData.
tEts.end()) {
259 const bool diagonal_block =
260 (row_type == col_type) && (row_side == col_side);
263 locK.resize(row_nb_dofs, col_nb_dofs,
false);
266 const int row_nb_gauss_pts = row_data.
getN().size1();
267 const int row_nb_base_functions = row_data.
getN().size2();
274 double density = massData.
rho0;
277 auto t_w = getFTensor0IntegrationWeight();
280 for (
int gg = 0; gg != row_nb_gauss_pts; gg++) {
285 double w = getVolume() * t_w;
287 for (
int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
289 for (
int col_bb = 0; col_bb != col_nb_dofs / 3; col_bb++) {
290 auto t_assemble = get_tensor2(locK, row_bb, col_bb);
291 t_assemble(
i,
j) += density * t_row_base_func * t_col_base_func * w;
306 if (row_type != col_type || row_side != col_side) {
307 translocK.resize(col_nb_dofs, row_nb_dofs,
false);
308 noalias(translocK) = trans(locK);
320 boost::shared_ptr<map<int, BlockData>>
331 const std::string row_field,
const std::string col_field,
332 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
333 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
334 boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
const double rho_n,
342 OpAssemble(
const std::string row_field,
const std::string col_field,
343 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
344 const char type,
bool symm =
false);
403 OpRhs_dx(
const std::string row_field,
const std::string col_field,
404 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
412 OpLhs_dx_dx(
const std::string row_field,
const std::string col_field,
413 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
427 OpAleRhs_dx(
const std::string row_field,
const std::string col_field,
428 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
437 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
452 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
472 const std::string row_field,
const std::string col_field,
473 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
474 boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
475 boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts,
476 const double rho_n,
const double rho_0);
496 const std::string row_field,
const std::string col_field,
497 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
498 boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
499 boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts,
500 const double rho_n,
const double rho_0);
514 OpAleRhs_dX(
const std::string row_field,
const std::string col_field,
515 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
524 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
539 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
550 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
551 :
OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false) {}
565 typedef boost::function<
577 const std::string row_field,
578 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
588 typedef boost::function<
600 const std::string row_field,
601 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
603 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr);
613 typedef boost::function<
625 const std::string row_field,
626 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
628 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr);
636 template <
class ELEMENT>
647 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
648 map<int, BlockData> &block_sets_ptr,
649 moab::Interface &post_proc_mesh,
650 std::vector<EntityHandle> &map_gauss_pts,
651 bool is_ale =
false,
bool is_field_disp =
true);
659 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr);
663 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
664 const std::string element_name,
const std::string x_field,
665 const std::string X_field,
const bool ale);
669 boost::shared_ptr<ForcesAndSourcesCore> fe_rhs_ptr,
670 boost::shared_ptr<map<int, BlockData>> block_sets_ptr,
671 const std::string x_field,
const std::string X_field,
672 const bool ale,
const bool field_disp,
674 boost::shared_ptr<DataAtIntegrationPts> data_at_pts =
nullptr);
678 const std::string x_field,
const std::string X_field,
679 const bool ale,
const bool field_disp,
687HookeElement::OpCalculateStrain<D>::OpCalculateStrain(
688 const std::string row_field,
const std::string col_field,
689 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
691 dataAtPts(data_at_pts) {
692 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE],
false);
696MoFEMErrorCode HookeElement::OpCalculateStrain<D>::doWork(
int row_side,
703 const int nb_integration_pts = getGaussPts().size2();
704 dataAtPts->smallStrainMat->resize(6, nb_integration_pts,
false);
705 auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
706 auto t_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hMat));
708 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
709 t_strain(
i,
j) = (t_h(
i,
j) || t_h(
j,
i)) / 2.;
725HookeElement::OpAleLhs_dx_dx<S>::OpAleLhs_dx_dx(
726 const std::string row_field,
const std::string col_field,
727 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
728 :
OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
739 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 1,
c + 0),
740 &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 2,
c + 0), &
m(r + 2,
c + 1),
750 double vol = getVolume();
753 auto t_w = getFTensor0IntegrationWeight();
757 const int row_nb_base_fun = row_data.
getN().size2();
764 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
765 auto &det_H = *dataAtPts->detHVec;
768 for (
int gg = 0; gg != nbIntegrationPts; ++gg) {
771 double a = t_w * vol * det_H[gg];
775 for (; rr != nbRows / 3; ++rr) {
778 auto t_m = get_tensor2(
K, 3 * rr, 0);
781 t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
787 t_rowD(
l,
j,
k) = t_D(
i,
j,
k,
l) * (
a * t_row_diff_base_pulled(
i));
793 for (
int cc = 0; cc != nbCols / 3; ++cc) {
796 t_col_diff_base_pulled(
j) = t_col_diff_base(
i) * t_invH(
i,
j);
799 t_m(
i,
j) += t_rowD(
i,
j,
k) * t_col_diff_base_pulled(
k);
812 for (; rr != row_nb_base_fun; ++rr)
825HookeElement::OpCalculateHomogeneousStiffness<S>::
826 OpCalculateHomogeneousStiffness(
827 const std::string row_field,
const std::string col_field,
828 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
829 boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
831 blockSetsPtr(block_sets_ptr), dataAtPts(data_at_pts) {
832 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE],
false);
836MoFEMErrorCode HookeElement::OpCalculateHomogeneousStiffness<S>::doWork(
840 for (
auto &
m : (*blockSetsPtr)) {
841 if (
m.second.tEts.find(getFEEntityHandle()) !=
m.second.tEts.end()) {
843 dataAtPts->stiffnessMat->resize(36, 1,
false);
846 const double young =
m.second.E;
847 const double poisson =
m.second.PoissonRatio;
850 const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
857 t_D(
i,
j,
k,
l) = 0.;
859 t_D(0, 0, 0, 0) = 1 - poisson;
860 t_D(1, 1, 1, 1) = 1 - poisson;
861 t_D(2, 2, 2, 2) = 1 - poisson;
863 t_D(0, 1, 0, 1) = 0.5 * (1 - 2 * poisson);
864 t_D(0, 2, 0, 2) = 0.5 * (1 - 2 * poisson);
865 t_D(1, 2, 1, 2) = 0.5 * (1 - 2 * poisson);
867 t_D(0, 0, 1, 1) = poisson;
868 t_D(1, 1, 0, 0) = poisson;
869 t_D(0, 0, 2, 2) = poisson;
870 t_D(2, 2, 0, 0) = poisson;
871 t_D(1, 1, 2, 2) = poisson;
872 t_D(2, 2, 1, 1) = poisson;
874 t_D(
i,
j,
k,
l) *= coefficient;
884HookeElement::OpCalculateStress<S>::OpCalculateStress(
885 const std::string row_field,
const std::string col_field,
886 boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
888 dataAtPts(data_at_pts) {
889 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE],
false);
893MoFEMErrorCode HookeElement::OpCalculateStress<S>::doWork(
int row_side,
898 const int nb_integration_pts = getGaussPts().size2();
899 auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
900 dataAtPts->cauchyStressMat->resize(6, nb_integration_pts,
false);
901 auto t_cauchy_stress =
902 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
913 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
914 t_cauchy_stress(
i,
j) = t_D(
i,
j,
k,
l) * t_strain(
k,
l);
923HookeElement::OpLhs_dx_dx<S>::OpLhs_dx_dx(
924 const std::string row_field,
const std::string col_field,
925 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
926 :
OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
937 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 1,
c + 0),
938 &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 2,
c + 0), &
m(r + 2,
c + 1),
948 double vol = getVolume();
951 auto t_w = getFTensor0IntegrationWeight();
955 const int row_nb_base_fun = row_data.
getN().size2();
963 for (
int gg = 0; gg != nbIntegrationPts; ++gg) {
966 double a = t_w * vol;
970 for (; rr != nbRows / 3; ++rr) {
973 auto t_m = get_tensor2(
K, 3 * rr, 0);
982 t_rowD(
l,
j,
k) = t_D(
i,
j,
k,
l) * (
a * t_row_diff_base(
i));
985 for (
int cc = 0; cc != nbCols / 3; ++cc) {
988 t_m(
i,
j) += t_rowD(
i,
j,
k) * t_col_diff_base(
k);
1001 for (; rr != row_nb_base_fun; ++rr)
1013HookeElement::OpAleLhs_dx_dX<S>::OpAleLhs_dx_dX(
1014 const std::string row_field,
const std::string col_field,
1015 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1016 :
OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false) {}
1027 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 1,
c + 0),
1028 &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 2,
c + 0), &
m(r + 2,
c + 1),
1040 double vol = getVolume();
1043 auto t_w = getFTensor0IntegrationWeight();
1047 const int row_nb_base_fun = row_data.
getN().size2();
1054 auto t_cauchy_stress =
1055 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1056 auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1057 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1058 auto &det_H = *dataAtPts->detHVec;
1061 for (
int gg = 0; gg != nbIntegrationPts; ++gg) {
1064 double a = t_w * vol * det_H[gg];
1067 t_F_dX(
i,
j,
k,
l) = -(t_h(
i,
m) * t_invH(
m,
k)) * t_invH(
l,
j);
1071 for (; rr != nbRows / 3; ++rr) {
1074 auto t_m = get_tensor2(
K, 3 * rr, 0);
1077 t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
1080 t_row_stress(
i) =
a * t_row_diff_base_pulled(
j) * t_cauchy_stress(
i,
j);
1083 t_row_diff_base_pulled_dX(
j,
k,
l) =
1084 -(t_invH(
i,
k) * t_row_diff_base(
i)) * t_invH(
l,
j);
1087 t_row_dX_stress(
i,
k,
l) =
1088 a * (t_row_diff_base_pulled_dX(
j,
k,
l) * t_cauchy_stress(
j,
i));
1091 t_row_D(
l,
j,
k) = (
a * t_row_diff_base_pulled(
i)) * t_D(
i,
j,
k,
l);
1096 t_row_stress_dX(
i,
j,
k) = 0;
1097 for (
int ii = 0; ii != 3; ++ii)
1098 for (
int mm = 0; mm != 3; ++mm)
1099 for (
int nn = 0; nn != 3; ++nn) {
1100 auto &
v = t_row_stress_dX(ii, mm, nn);
1101 for (
int kk = 0; kk != 3; ++kk)
1102 for (
int ll = 0; ll != 3; ++ll)
1103 v += t_row_D(ii, kk, ll) * t_F_dX(kk, ll, mm, nn);
1110 for (
int cc = 0; cc != nbCols / 3; ++cc) {
1112 t_m(
i,
k) += t_row_stress(
i) * (t_invH(
j,
k) * t_col_diff_base(
j));
1113 t_m(
i,
k) += t_row_dX_stress(
i,
k,
l) * t_col_diff_base(
l);
1114 t_m(
i,
k) += t_row_stress_dX(
i,
k,
l) * t_col_diff_base(
l);
1127 for (; rr != row_nb_base_fun; ++rr)
1142HookeElement::OpAleLhs_dX_dX<S>::OpAleLhs_dX_dX(
1143 const std::string row_field,
const std::string col_field,
1144 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1145 :
OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
1156 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 1,
c + 0),
1157 &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 2,
c + 0), &
m(r + 2,
c + 1),
1169 double vol = getVolume();
1172 auto t_w = getFTensor0IntegrationWeight();
1176 const int row_nb_base_fun = row_data.
getN().size2();
1182 auto t_cauchy_stress =
1183 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1184 auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
1185 auto t_eshelby_stress =
1186 getFTensor2FromMat<3, 3>(*dataAtPts->eshelbyStressMat);
1187 auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1188 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1189 auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
1190 auto &det_H = *dataAtPts->detHVec;
1193 for (
int gg = 0; gg != nbIntegrationPts; ++gg) {
1196 double a = t_w * vol * det_H[gg];
1199 t_F_dX(
i,
j,
k,
l) = -(t_h(
i,
m) * t_invH(
m,
k)) * t_invH(
l,
j);
1202 t_D_strain_dX(
i,
j,
m,
n) = 0.;
1203 for (
int ii = 0; ii != 3; ++ii)
1204 for (
int jj = 0; jj != 3; ++jj)
1205 for (
int ll = 0; ll != 3; ++ll)
1206 for (
int kk = 0; kk != 3; ++kk) {
1207 auto &
v = t_D_strain_dX(ii, jj, kk, ll);
1208 for (
int mm = 0; mm != 3; ++mm)
1209 for (
int nn = 0; nn != 3; ++nn)
1210 v += t_D(ii, jj, mm, nn) * t_F_dX(mm, nn, kk, ll);
1214 t_eshelby_stress_dX(
i,
j,
m,
n) = t_F(
k,
i) * t_D_strain_dX(
k,
j,
m,
n);
1216 for (
int ii = 0; ii != 3; ++ii)
1217 for (
int jj = 0; jj != 3; ++jj)
1218 for (
int mm = 0; mm != 3; ++mm)
1219 for (
int nn = 0; nn != 3; ++nn) {
1220 auto &
v = t_eshelby_stress_dX(ii, jj, mm, nn);
1221 for (
int kk = 0; kk != 3; ++kk)
1222 v += t_F_dX(kk, ii, mm, nn) * t_cauchy_stress(kk, jj);
1225 t_eshelby_stress_dX(
i,
j,
k,
l) *= -1;
1228 t_energy_dX(
k,
l) = t_F_dX(
i,
j,
k,
l) * t_cauchy_stress(
i,
j);
1229 t_energy_dX(
k,
l) +=
1230 (t_strain(
m,
n) * t_D(
m,
n,
i,
j)) * t_F_dX(
i,
j,
k,
l);
1231 t_energy_dX(
k,
l) /= 2.;
1233 for (
int kk = 0; kk != 3; ++kk)
1234 for (
int ll = 0; ll != 3; ++ll) {
1235 auto v = t_energy_dX(kk, ll);
1236 for (
int ii = 0; ii != 3; ++ii)
1237 t_eshelby_stress_dX(ii, ii, kk, ll) +=
v;
1242 for (; rr != nbRows / 3; ++rr) {
1245 auto t_m = get_tensor2(
K, 3 * rr, 0);
1248 t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
1251 t_row_stress(
i) =
a * t_row_diff_base_pulled(
j) * t_eshelby_stress(
i,
j);
1254 t_row_diff_base_pulled_dX(
j,
k,
l) =
1255 -(t_row_diff_base(
i) * t_invH(
i,
k)) * t_invH(
l,
j);
1258 t_row_dX_stress(
i,
k,
l) =
1259 a * (t_row_diff_base_pulled_dX(
j,
k,
l) * t_eshelby_stress(
i,
j));
1262 t_row_stress_dX(
i,
m,
n) =
1263 a * t_row_diff_base_pulled(
j) * t_eshelby_stress_dX(
i,
j,
m,
n);
1269 for (
int cc = 0; cc != nbCols / 3; ++cc) {
1271 t_m(
i,
k) += t_row_stress(
i) * (t_invH(
j,
k) * t_col_diff_base(
j));
1272 t_m(
i,
k) += t_row_dX_stress(
i,
k,
l) * t_col_diff_base(
l);
1273 t_m(
i,
k) += t_row_stress_dX(
i,
k,
l) * t_col_diff_base(
l);
1286 for (; rr != row_nb_base_fun; ++rr)
1304HookeElement::OpAleLhsPre_dX_dx<S>::OpAleLhsPre_dX_dx(
1305 const std::string row_field,
const std::string col_field,
1306 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1308 dataAtPts(data_at_pts) {
1309 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE],
false);
1313MoFEMErrorCode HookeElement::OpAleLhsPre_dX_dx<S>::doWork(
int row_side,
1318 const int nb_integration_pts = row_data.
getN().size1();
1320 auto get_eshelby_stress_dx = [
this, nb_integration_pts]() {
1322 t_eshelby_stress_dx;
1323 dataAtPts->eshelbyStress_dx->resize(81, nb_integration_pts,
false);
1325 for (
int ii = 0; ii != 3; ++ii)
1326 for (
int jj = 0; jj != 3; ++jj)
1327 for (
int kk = 0; kk != 3; ++kk)
1328 for (
int ll = 0; ll != 3; ++ll)
1329 t_eshelby_stress_dx.ptr(ii, jj, kk, ll) =
1330 &(*dataAtPts->eshelbyStress_dx)(mm++, 0);
1331 return t_eshelby_stress_dx;
1334 auto t_eshelby_stress_dx = get_eshelby_stress_dx();
1347 auto t_cauchy_stress =
1348 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1349 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1350 auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
1352 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1354 t_eshelby_stress_dx(
i,
j,
m,
n) =
1355 (t_F(
k,
i) * t_D(
k,
j,
m,
l)) * t_invH(
n,
l);
1356 for (
int ii = 0; ii != 3; ++ii)
1357 for (
int jj = 0; jj != 3; ++jj)
1358 for (
int mm = 0; mm != 3; ++mm)
1359 for (
int nn = 0; nn != 3; ++nn) {
1360 auto &
v = t_eshelby_stress_dx(ii, jj, mm, nn);
1361 v += t_invH(nn, ii) * t_cauchy_stress(mm, jj);
1363 t_eshelby_stress_dx(
i,
j,
k,
l) *= -1;
1366 t_energy_dx(
m,
n) = t_invH(
n,
j) * t_cauchy_stress(
m,
j);
1368 for (
int mm = 0; mm != 3; ++mm)
1369 for (
int nn = 0; nn != 3; ++nn) {
1370 auto v = t_energy_dx(mm, nn);
1371 for (
int ii = 0; ii != 3; ++ii)
1372 t_eshelby_stress_dx(ii, ii, mm, nn) +=
v;
1378 ++t_eshelby_stress_dx;
1385template <
class ELEMENT>
1386HookeElement::OpPostProcHookeElement<ELEMENT>::OpPostProcHookeElement(
1387 const string row_field, boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1388 map<int, BlockData> &block_sets_ptr, moab::Interface &post_proc_mesh,
1389 std::vector<EntityHandle> &map_gauss_pts,
bool is_ale,
bool is_field_disp)
1391 dataAtPts(data_at_pts), blockSetsPtr(block_sets_ptr),
1392 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts), isALE(is_ale),
1393 isFieldDisp(is_field_disp) {}
1395template <
class ELEMENT>
1396MoFEMErrorCode HookeElement::OpPostProcHookeElement<ELEMENT>::doWork(
1400 if (type != MBVERTEX) {
1404 auto tensor_to_tensor = [](
const auto &t1,
auto &t2) {
1405 t2(0, 0) = t1(0, 0);
1406 t2(1, 1) = t1(1, 1);
1407 t2(2, 2) = t1(2, 2);
1408 t2(0, 1) = t2(1, 0) = t1(1, 0);
1409 t2(0, 2) = t2(2, 0) = t1(2, 0);
1410 t2(1, 2) = t2(2, 1) = t1(2, 1);
1413 std::array<double, 9> def_val;
1416 auto make_tag = [&](
auto name,
auto size) {
1418 CHKERR postProcMesh.tag_get_handle(name, size, MB_TYPE_DOUBLE,
th,
1419 MB_TAG_CREAT | MB_TAG_SPARSE,
1424 auto th_stress = make_tag(
"STRESS", 9);
1425 auto th_psi = make_tag(
"ENERGY", 1);
1427 const int nb_integration_pts = mapGaussPts.size();
1434 auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1435 auto t_H = getFTensor2FromMat<3, 3>(*dataAtPts->HMat);
1437 dataAtPts->stiffnessMat->resize(36, 1,
false);
1444 if (type == MBTRI || type == MBQUAD) {
1446 auto &m_field = this->getPtrFE()->mField;
1447 CHKERR m_field.get_moab().get_adjacencies(&ent, 1, 3,
false, ents,
1448 moab::Interface::UNION);
1452 "Could not find a 3D element adjacent to a given face element");
1454 ent_3d = ents.front();
1457 bool found_block =
false;
1458 for (
auto &
m : (blockSetsPtr)) {
1459 if (
m.second.tEts.find(ent_3d) !=
m.second.tEts.end()) {
1460 const double young =
m.second.E;
1461 const double poisson =
m.second.PoissonRatio;
1462 const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
1464 t_D(
i,
j,
k,
l) = 0.;
1465 t_D(0, 0, 0, 0) = t_D(1, 1, 1, 1) = t_D(2, 2, 2, 2) = 1 - poisson;
1466 t_D(0, 1, 0, 1) = t_D(0, 2, 0, 2) = t_D(1, 2, 1, 2) =
1467 0.5 * (1 - 2 * poisson);
1468 t_D(0, 0, 1, 1) = t_D(1, 1, 0, 0) = t_D(0, 0, 2, 2) = t_D(2, 2, 0, 0) =
1469 t_D(1, 1, 2, 2) = t_D(2, 2, 1, 1) = poisson;
1470 t_D(
i,
j,
k,
l) *= coefficient;
1478 "Element not found in any of material blocksets");
1488 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1497 t_small_strain_symm(
i,
j) = (t_h(
i,
j) || t_h(
j,
i)) / 2.;
1501 t_F(
i,
j) = t_h(
i,
k) * t_invH(
k,
j);
1502 t_small_strain_symm(
i,
j) = (t_F(
i,
j) || t_F(
j,
i)) / 2.;
1506 t_small_strain_symm(0, 0) -= 1;
1507 t_small_strain_symm(1, 1) -= 1;
1508 t_small_strain_symm(2, 2) -= 1;
1511 t_stress_symm(
i,
j) = t_D(
i,
j,
k,
l) * t_small_strain_symm(
k,
l);
1512 tensor_to_tensor(t_stress_symm, t_stress);
1514 const double psi = 0.5 * t_stress_symm(
i,
j) * t_small_strain_symm(
i,
j);
1516 CHKERR postProcMesh.tag_set_data(th_psi, &mapGaussPts[gg], 1, &psi);
1517 CHKERR postProcMesh.tag_set_data(th_stress, &mapGaussPts[gg], 1,
1527HookeElement::OpAnalyticalInternalStrain_dx<S>::OpAnalyticalInternalStrain_dx(
1528 const std::string row_field,
1529 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1530 StrainFunction strain_fun)
1531 :
OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1532 strainFun(strain_fun) {}
1536HookeElement::OpAnalyticalInternalStrain_dx<S>::iNtegrate(
EntData &row_data) {
1545 &
v(r + 0), &
v(r + 1), &
v(r + 2));
1548 const int nb_integration_pts = getGaussPts().size2();
1549 auto t_coords = getFTensor1CoordsAtGaussPts();
1552 double vol = getVolume();
1553 auto t_w = getFTensor0IntegrationWeight();
1555 nF.resize(nbRows,
false);
1565 const int row_nb_base_fun = row_data.
getN().size2();
1567 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
1569 auto t_fun_strain = strainFun(t_coords);
1571 t_stress(
i,
j) = -t_D(
i,
j,
k,
l) * t_fun_strain(
k,
l);
1574 double a = t_w * vol;
1576 auto t_nf = get_tensor1(nF, 0);
1579 for (; rr != nbRows / 3; ++rr) {
1580 t_nf(
i) +=
a * t_row_diff_base(
j) * t_stress(
i,
j);
1585 for (; rr != row_nb_base_fun; ++rr)
1597HookeElement::OpAnalyticalInternalAleStrain_dX<S>::
1598 OpAnalyticalInternalAleStrain_dX(
1599 const std::string row_field,
1600 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
1601 StrainFunction strain_fun,
1602 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr)
1603 :
OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1604 strainFun(strain_fun), matPosAtPtsPtr(mat_pos_at_pts_ptr) {}
1607MoFEMErrorCode HookeElement::OpAnalyticalInternalAleStrain_dX<S>::iNtegrate(
1617 &
v(r + 0), &
v(r + 1), &
v(r + 2));
1620 const int nb_integration_pts = getGaussPts().size2();
1622 auto get_coords = [&]() {
return getFTensor1FromMat<3>(*matPosAtPtsPtr); };
1623 auto t_coords = get_coords();
1626 double vol = getVolume();
1627 auto t_w = getFTensor0IntegrationWeight();
1629 nF.resize(nbRows,
false);
1636 auto t_F = getFTensor2FromMat<3, 3>(*(dataAtPts->FMat));
1637 auto &det_H = *dataAtPts->detHVec;
1638 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1642 const int row_nb_base_fun = row_data.
getN().size2();
1644 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
1646 auto t_fun_strain = strainFun(t_coords);
1648 t_stress(
i,
j) = -t_D(
i,
j,
k,
l) * t_fun_strain(
k,
l);
1650 t_eshelby_stress(
i,
j) = -t_F(
k,
i) * t_stress(
k,
j);
1653 double a = t_w * vol * det_H[gg];
1655 auto t_nf = get_tensor1(nF, 0);
1658 for (; rr != nbRows / 3; ++rr) {
1660 t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
1661 t_nf(
i) +=
a * t_row_diff_base_pulled(
j) * t_eshelby_stress(
i,
j);
1666 for (; rr != row_nb_base_fun; ++rr)
1680HookeElement::OpAnalyticalInternalAleStrain_dx<S>::
1681 OpAnalyticalInternalAleStrain_dx(
1682 const std::string row_field,
1683 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
1684 StrainFunction strain_fun,
1685 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr)
1686 :
OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1687 strainFun(strain_fun), matPosAtPtsPtr(mat_pos_at_pts_ptr) {}
1690MoFEMErrorCode HookeElement::OpAnalyticalInternalAleStrain_dx<S>::iNtegrate(
1700 &
v(r + 0), &
v(r + 1), &
v(r + 2));
1703 const int nb_integration_pts = getGaussPts().size2();
1705 auto get_coords = [&]() {
return getFTensor1FromMat<3>(*matPosAtPtsPtr); };
1706 auto t_coords = get_coords();
1709 double vol = getVolume();
1710 auto t_w = getFTensor0IntegrationWeight();
1712 nF.resize(nbRows,
false);
1719 auto &det_H = *dataAtPts->detHVec;
1720 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1724 const int row_nb_base_fun = row_data.
getN().size2();
1726 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
1728 auto t_fun_strain = strainFun(t_coords);
1730 t_stress(
i,
j) = -t_D(
i,
j,
k,
l) * t_fun_strain(
k,
l);
1733 double a = t_w * vol * det_H[gg];
1735 auto t_nf = get_tensor1(nF, 0);
1738 for (; rr != nbRows / 3; ++rr) {
1740 t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
1741 t_nf(
i) +=
a * t_row_diff_base_pulled(
j) * t_stress(
i,
j);
1746 for (; rr != row_nb_base_fun; ++rr)
static MoFEMErrorCode setOperators(boost::shared_ptr< ForcesAndSourcesCore > fe_lhs_ptr, boost::shared_ptr< ForcesAndSourcesCore > fe_rhs_ptr, boost::shared_ptr< map< int, BlockData > > block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, const EntityType type=MBTET, boost::shared_ptr< DataAtIntegrationPts > data_at_pts=nullptr)
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData > > &block_sets_ptr)
static MoFEMErrorCode addElasticElement(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData > > &block_sets_ptr, const std::string element_name, const std::string x_field, const std::string X_field, const bool ale)
static MoFEMErrorCode calculateEnergy(DM dm, boost::shared_ptr< map< int, BlockData > > block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, SmartPetscObj< Vec > &v_energy_ptr)
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#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_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
UBlasVector< double > VectorDouble
UBlasVector< int > VectorInt
auto type_from_handle(const EntityHandle h)
get type from entity handle
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 determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
data for calculation inertia forces
Range tEts
elements in block set
double rho0
reference density
VectorDouble a0
constant acceleration
structure grouping operators and data used for calculation of mass (convective) element \ nonlinear_e...
boost::shared_ptr< MatrixDouble > HMat
boost::shared_ptr< MatrixDouble > invHMat
Range forcesOnlyOnEntitiesRow
boost::shared_ptr< MatrixDouble > smallStrainMat
boost::shared_ptr< MatrixDouble > cauchyStressMat
boost::shared_ptr< MatrixDouble > eshelbyStress_dx
boost::shared_ptr< MatrixDouble > FMat
Range forcesOnlyOnEntitiesCol
boost::shared_ptr< MatrixDouble > eshelbyStressMat
boost::shared_ptr< MatrixDouble > stiffnessMat
boost::shared_ptr< MatrixDouble > hMat
boost::shared_ptr< VectorDouble > energyVec
boost::shared_ptr< VectorDouble > detHVec
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.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of dofs on entity.
intrusive_ptr for managing petsc objects
Volume finite element base.
data for calculation heat conductivity and heat capacity elements
Range forcesOnlyOnEntitiesRow
Range forcesOnlyOnEntitiesCol
Range tEts
constrains elements in block set
structure grouping operators and data used for calculation of nonlinear elastic element
OpAleLhs_dX_dX(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate tangent stiffness for material momentum.
OpAleLhs_dX_dx(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate tangent stiffness for material momentum.
OpAleLhs_dx_dX(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate tangent stiffness for spatial momentum.
OpAleLhs_dx_dx(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate B^T D B operator.
OpAleLhsPre_dX_dx(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate tangent stiffness for material momentum.
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
OpAleLhsWithDensity_dX_dX(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts, boost::shared_ptr< VectorDouble > rho_at_gauss_pts, boost::shared_ptr< MatrixDouble > rho_grad_at_gauss_pts, const double rho_n, const double rho_0)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate tangent stiffness for spatial momentum.
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
OpAleLhsWithDensity_dx_dX(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts, boost::shared_ptr< VectorDouble > rho_at_gauss_pts, boost::shared_ptr< MatrixDouble > rho_grad_at_gauss_pts, const double rho_n, const double rho_0)
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
MoFEMErrorCode iNtegrate(EntData &row_data)
OpAleRhs_dX(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
OpAleRhs_dx(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
MoFEMErrorCode iNtegrate(EntData &row_data)
MoFEMErrorCode iNtegrate(EntData &row_data)
OpAnalyticalInternalAleStrain_dX(const std::string row_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts, StrainFunction strain_fun, boost::shared_ptr< MatrixDouble > mat_pos_at_pts_ptr)
boost::shared_ptr< MatrixDouble > matPosAtPtsPtr
boost::function< FTensor::Tensor2_symmetric< double, 3 >(FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &t_coords) > StrainFunction
OpAnalyticalInternalAleStrain_dx(const std::string row_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts, StrainFunction strain_fun, boost::shared_ptr< MatrixDouble > mat_pos_at_pts_ptr)
boost::shared_ptr< MatrixDouble > matPosAtPtsPtr
MoFEMErrorCode iNtegrate(EntData &row_data)
boost::function< FTensor::Tensor2_symmetric< double, 3 >(FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &t_coords) > StrainFunction
OpAnalyticalInternalStrain_dx(const std::string row_field, boost::shared_ptr< DataAtIntegrationPts > data_at_pts, StrainFunction strain_fun)
boost::function< FTensor::Tensor2_symmetric< double, 3 >(FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, 3 > &t_coords) > StrainFunction
MoFEMErrorCode iNtegrate(EntData &row_data)
int nbCols
number if dof on column
bool isDiag
true if this block is on diagonal
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
virtual MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
Assemble local entity block matrix.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Do calculations for give operator.
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
int nbIntegrationPts
number of integration points
int nbRows
number of dofs on rows
virtual MoFEMErrorCode iNtegrate(EntData &row_data)
MoFEMErrorCode aSsemble(EntData &row_data)
Assemble local entity right-hand vector.
OpAssemble(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_at_pts, const char type, bool symm=false)
OpCalculateEnergy(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_at_pts, SmartPetscObj< Vec > ghost_vec=SmartPetscObj< Vec >())
SmartPetscObj< Vec > ghostVec
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
OpCalculateEshelbyStress(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_at_pts)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< map< int, BlockData > > blockSetsPtr
OpCalculateHomogeneousStiffness(const std::string row_field, const std::string col_field, boost::shared_ptr< map< int, BlockData > > &block_sets_ptr, boost::shared_ptr< DataAtIntegrationPts > data_at_pts)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Assemble mass matrix for elastic element TODO: CHANGE FORMULA *.
OpCalculateMassMatrix(const std::string row_field, const std::string col_field, BlockData &data, MassBlockData &mass_data, boost::shared_ptr< DataAtIntegrationPts > &common_data, bool symm=true)
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
boost::shared_ptr< DataAtIntegrationPts > commonData
const double rHo0
p_0 reference density in E(p) = E * (p / p_0)^n
boost::shared_ptr< map< int, BlockData > > blockSetsPtr
OpCalculateStiffnessScaledByDensityField(const std::string row_field, const std::string col_field, boost::shared_ptr< map< int, BlockData > > &block_sets_ptr, boost::shared_ptr< DataAtIntegrationPts > data_at_pts, boost::shared_ptr< VectorDouble > rho_at_gauss_pts, const double rho_n, const double rho_0)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
const double rhoN
exponent n in E(p) = E * (p / p_0)^n
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateStrainAle(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
OpCalculateStrain(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateStress(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_at_pts)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
OpLhs_dx_dx(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate B^T D B operator.
std::vector< EntityHandle > & mapGaussPts
map< int, BlockData > & blockSetsPtr
moab::Interface & postProcMesh
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
OpPostProcHookeElement(const string row_field, boost::shared_ptr< DataAtIntegrationPts > data_at_pts, map< int, BlockData > &block_sets_ptr, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, bool is_ale=false, bool is_field_disp=true)
OpRhs_dx(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
MoFEMErrorCode iNtegrate(EntData &row_data)