31 #ifndef __HOOKE_ELEMENT_HPP
32 #define __HOOKE_ELEMENT_HPP
34 #ifndef __BASICFINITEELEMENTS_HPP__
38 #ifndef __NONLINEAR_ELASTIC_HPP
69 #ifndef __CONVECTIVE_MASS_ELEMENT_HPP
95 boost::shared_ptr<MatrixDouble>
hMat;
96 boost::shared_ptr<MatrixDouble>
FMat;
98 boost::shared_ptr<MatrixDouble>
HMat;
111 smallStrainMat = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
112 hMat = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
113 FMat = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
115 HMat = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
116 detHVec = boost::shared_ptr<VectorDouble>(
new VectorDouble());
117 invHMat = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
119 cauchyStressMat = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
120 stiffnessMat = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
121 energyVec = boost::shared_ptr<VectorDouble>(
new VectorDouble());
122 eshelbyStressMat = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
123 stiffnessMat = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
125 eshelbyStress_dx = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
132 template <
bool D = false>
136 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
147 const std::string col_field,
148 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
156 #define MAT_TO_DDG(SM) \
157 &(*SM)(0, 0), &(*SM)(1, 0), &(*SM)(2, 0), &(*SM)(3, 0), &(*SM)(4, 0), \
158 &(*SM)(5, 0), &(*SM)(6, 0), &(*SM)(7, 0), &(*SM)(8, 0), &(*SM)(9, 0), \
159 &(*SM)(10, 0), &(*SM)(11, 0), &(*SM)(12, 0), &(*SM)(13, 0), \
160 &(*SM)(14, 0), &(*SM)(15, 0), &(*SM)(16, 0), &(*SM)(17, 0), \
161 &(*SM)(18, 0), &(*SM)(19, 0), &(*SM)(20, 0), &(*SM)(21, 0), \
162 &(*SM)(22, 0), &(*SM)(23, 0), &(*SM)(24, 0), &(*SM)(25, 0), \
163 &(*SM)(26, 0), &(*SM)(27, 0), &(*SM)(28, 0), &(*SM)(29, 0), \
164 &(*SM)(30, 0), &(*SM)(31, 0), &(*SM)(32, 0), &(*SM)(33, 0), \
165 &(*SM)(34, 0), &(*SM)(35, 0)
170 boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
181 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
182 Vec ghost_vec = PETSC_NULL);
196 const std::string row_field,
const std::string col_field,
197 boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
209 const std::string row_field,
const std::string col_field,
210 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
211 boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
216 boost::shared_ptr<map<int, BlockData>>
239 const std::string col_field,
BlockData &data,
241 boost::shared_ptr<DataAtIntegrationPts> &common_data,
244 row_field, col_field, OPROWCOL, symm),
245 commonData(common_data), dAta(data), massData(mass_data) {}
247 PetscErrorCode
doWork(
int row_side,
int col_side, EntityType row_type,
255 &
m(3 *
r + 0, 3 *
c + 0), &
m(3 *
r + 0, 3 *
c + 1),
256 &
m(3 *
r + 0, 3 *
c + 2), &
m(3 *
r + 1, 3 *
c + 0),
257 &
m(3 *
r + 1, 3 *
c + 1), &
m(3 *
r + 1, 3 *
c + 2),
258 &
m(3 *
r + 2, 3 *
c + 0), &
m(3 *
r + 2, 3 *
c + 1),
259 &
m(3 *
r + 2, 3 *
c + 2));
262 const int row_nb_dofs = row_data.
getIndices().size();
265 const int col_nb_dofs = col_data.
getIndices().size();
268 if (dAta.tEts.find(getFEEntityHandle()) == dAta.tEts.end()) {
271 if (massData.
tEts.find(getFEEntityHandle()) == massData.
tEts.end()) {
275 const bool diagonal_block =
276 (row_type == col_type) && (row_side == col_side);
279 locK.resize(row_nb_dofs, col_nb_dofs,
false);
282 const int row_nb_gauss_pts = row_data.
getN().size1();
283 const int row_nb_base_functions = row_data.
getN().size2();
290 double density = massData.
rho0;
293 auto t_w = getFTensor0IntegrationWeight();
296 for (
int gg = 0; gg != row_nb_gauss_pts; gg++) {
301 double w = getVolume() * t_w;
303 for (
int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
305 for (
int col_bb = 0; col_bb != col_nb_dofs / 3; col_bb++) {
306 auto t_assemble = get_tensor2(locK, row_bb, col_bb);
307 t_assemble(
i,
j) += density * t_row_base_func * t_col_base_func *
w;
322 if (row_type != col_type || row_side != col_side) {
323 translocK.resize(col_nb_dofs, row_nb_dofs,
false);
324 noalias(translocK) = trans(locK);
336 boost::shared_ptr<map<int, BlockData>>
347 const std::string row_field,
const std::string col_field,
348 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
349 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
350 boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
const double rho_n,
358 OpAssemble(
const std::string row_field,
const std::string col_field,
359 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
360 const char type,
bool symm =
false);
372 MoFEMErrorCode doWork(
int row_side,
int col_side, EntityType row_type,
373 EntityType col_type,
EntData &row_data,
419 OpRhs_dx(
const std::string row_field,
const std::string col_field,
420 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
428 OpLhs_dx_dx(
const std::string row_field,
const std::string col_field,
429 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
443 OpAleRhs_dx(
const std::string row_field,
const std::string col_field,
444 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
452 OpAleLhs_dx_dx(
const std::string row_field,
const std::string col_field,
453 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
467 OpAleLhs_dx_dX(
const std::string row_field,
const std::string col_field,
468 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
488 const std::string row_field,
const std::string col_field,
489 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
490 boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
491 boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts,
492 const double rho_n,
const double rho_0);
512 const std::string row_field,
const std::string col_field,
513 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
514 boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
515 boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts,
516 const double rho_n,
const double rho_0);
530 OpAleRhs_dX(
const std::string row_field,
const std::string col_field,
531 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
539 OpAleLhs_dX_dX(
const std::string row_field,
const std::string col_field,
540 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
555 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
566 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
567 :
OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false) {}
581 typedef boost::function<
593 const std::string row_field,
594 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
604 typedef boost::function<
616 const std::string row_field,
617 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
619 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr);
629 typedef boost::function<
641 const std::string row_field,
642 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
644 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr);
652 template <
class ELEMENT>
663 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
664 map<int, BlockData> &block_sets_ptr,
666 std::vector<EntityHandle> &map_gauss_pts,
667 bool is_ale =
false,
bool is_field_disp =
true);
675 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr);
679 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
680 const std::string element_name,
const std::string x_field,
681 const std::string X_field,
const bool ale);
685 boost::shared_ptr<ForcesAndSourcesCore> fe_rhs_ptr,
686 boost::shared_ptr<map<int, BlockData>> block_sets_ptr,
687 const std::string x_field,
const std::string X_field,
688 const bool ale,
const bool field_disp,
689 const EntityType
type = MBTET,
690 boost::shared_ptr<DataAtIntegrationPts> data_at_pts =
nullptr);
694 const std::string x_field,
const std::string X_field,
695 const bool ale,
const bool field_disp,
Vec *v_energy_ptr);
702 HookeElement::OpCalculateStrain<D>::OpCalculateStrain(
703 const std::string row_field,
const std::string col_field,
704 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
706 dataAtPts(data_at_pts) {
718 const int nb_integration_pts = getGaussPts().size2();
719 dataAtPts->smallStrainMat->resize(6, nb_integration_pts,
false);
720 auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
721 auto t_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hMat));
723 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
724 t_strain(
i,
j) = (t_h(
i,
j) || t_h(
j,
i)) / 2.;
740 HookeElement::OpAleLhs_dx_dx<S>::OpAleLhs_dx_dx(
741 const std::string row_field,
const std::string col_field,
742 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
743 :
OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
754 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2), &
m(
r + 1,
c + 0),
755 &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2), &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1),
765 double vol = getVolume();
768 auto t_w = getFTensor0IntegrationWeight();
772 const int row_nb_base_fun = row_data.
getN().size2();
779 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
780 auto &det_H = *dataAtPts->detHVec;
783 for (
int gg = 0; gg != nbIntegrationPts; ++gg) {
786 double a = t_w * vol * det_H[gg];
790 for (; rr != nbRows / 3; ++rr) {
793 auto t_m = get_tensor2(
K, 3 * rr, 0);
796 t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
802 t_rowD(
l,
j,
k) = t_D(
i,
j,
k,
l) * (a * t_row_diff_base_pulled(
i));
808 for (
int cc = 0; cc != nbCols / 3; ++cc) {
811 t_col_diff_base_pulled(
j) = t_col_diff_base(
i) * t_invH(
i,
j);
814 t_m(
i,
j) += t_rowD(
i,
j,
k) * t_col_diff_base_pulled(
k);
827 for (; rr != row_nb_base_fun; ++rr)
840 HookeElement::OpCalculateHomogeneousStiffness<S>::
841 OpCalculateHomogeneousStiffness(
842 const std::string row_field,
const std::string col_field,
843 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
844 boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
846 blockSetsPtr(block_sets_ptr), dataAtPts(data_at_pts) {
852 int row_side, EntityType row_type,
EntData &row_data) {
855 for (
auto &
m : (*blockSetsPtr)) {
856 if (
m.second.tEts.find(getFEEntityHandle()) !=
m.second.tEts.end()) {
858 dataAtPts->stiffnessMat->resize(36, 1,
false);
861 const double young =
m.second.E;
862 const double poisson =
m.second.PoissonRatio;
865 const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
872 t_D(
i,
j,
k,
l) = 0.;
874 t_D(0, 0, 0, 0) = 1 - poisson;
875 t_D(1, 1, 1, 1) = 1 - poisson;
876 t_D(2, 2, 2, 2) = 1 - poisson;
878 t_D(0, 1, 0, 1) = 0.5 * (1 - 2 * poisson);
879 t_D(0, 2, 0, 2) = 0.5 * (1 - 2 * poisson);
880 t_D(1, 2, 1, 2) = 0.5 * (1 - 2 * poisson);
882 t_D(0, 0, 1, 1) = poisson;
883 t_D(1, 1, 0, 0) = poisson;
884 t_D(0, 0, 2, 2) = poisson;
885 t_D(2, 2, 0, 0) = poisson;
886 t_D(1, 1, 2, 2) = poisson;
887 t_D(2, 2, 1, 1) = poisson;
889 t_D(
i,
j,
k,
l) *= coefficient;
899 HookeElement::OpCalculateStress<S>::OpCalculateStress(
900 const std::string row_field,
const std::string col_field,
901 boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
903 dataAtPts(data_at_pts) {
913 const int nb_integration_pts = getGaussPts().size2();
914 auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
915 dataAtPts->cauchyStressMat->resize(6, nb_integration_pts,
false);
916 auto t_cauchy_stress =
917 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
928 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
929 t_cauchy_stress(
i,
j) = t_D(
i,
j,
k,
l) * t_strain(
k,
l);
938 HookeElement::OpLhs_dx_dx<S>::OpLhs_dx_dx(
939 const std::string row_field,
const std::string col_field,
940 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
941 :
OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
952 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2), &
m(
r + 1,
c + 0),
953 &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2), &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1),
963 double vol = getVolume();
966 auto t_w = getFTensor0IntegrationWeight();
970 const int row_nb_base_fun = row_data.
getN().size2();
978 for (
int gg = 0; gg != nbIntegrationPts; ++gg) {
981 double a = t_w * vol;
982 if (getHoGaussPtsDetJac().size()) {
983 a *= getHoGaussPtsDetJac()[gg];
988 for (; rr != nbRows / 3; ++rr) {
991 auto t_m = get_tensor2(
K, 3 * rr, 0);
1000 t_rowD(
l,
j,
k) = t_D(
i,
j,
k,
l) * (a * t_row_diff_base(
i));
1003 for (
int cc = 0; cc != nbCols / 3; ++cc) {
1006 t_m(
i,
j) += t_rowD(
i,
j,
k) * t_col_diff_base(
k);
1019 for (; rr != row_nb_base_fun; ++rr)
1031 HookeElement::OpAleLhs_dx_dX<S>::OpAleLhs_dx_dX(
1032 const std::string row_field,
const std::string col_field,
1033 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1034 :
OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false) {}
1045 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2), &
m(
r + 1,
c + 0),
1046 &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2), &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1),
1058 double vol = getVolume();
1061 auto t_w = getFTensor0IntegrationWeight();
1065 const int row_nb_base_fun = row_data.
getN().size2();
1072 auto t_cauchy_stress =
1073 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1074 auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1075 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1076 auto &det_H = *dataAtPts->detHVec;
1079 for (
int gg = 0; gg != nbIntegrationPts; ++gg) {
1082 double a = t_w * vol * det_H[gg];
1085 t_F_dX(
i,
j,
k,
l) = -(t_h(
i,
m) * t_invH(
m,
k)) * t_invH(
l,
j);
1089 for (; rr != nbRows / 3; ++rr) {
1092 auto t_m = get_tensor2(
K, 3 * rr, 0);
1095 t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
1098 t_row_stress(
i) = a * t_row_diff_base_pulled(
j) * t_cauchy_stress(
i,
j);
1101 t_row_diff_base_pulled_dX(
j,
k,
l) =
1102 -(t_invH(
i,
k) * t_row_diff_base(
i)) * t_invH(
l,
j);
1105 t_row_dX_stress(
i,
k,
l) =
1106 a * (t_row_diff_base_pulled_dX(
j,
k,
l) * t_cauchy_stress(
j,
i));
1109 t_row_D(
l,
j,
k) = (a * t_row_diff_base_pulled(
i)) * t_D(
i,
j,
k,
l);
1114 t_row_stress_dX(
i,
j,
k) = 0;
1115 for (
int ii = 0; ii != 3; ++ii)
1116 for (
int mm = 0; mm != 3; ++mm)
1117 for (
int nn = 0; nn != 3; ++nn) {
1118 auto &v = t_row_stress_dX(ii, mm, nn);
1119 for (
int kk = 0; kk != 3; ++kk)
1120 for (
int ll = 0; ll != 3; ++ll)
1121 v += t_row_D(ii, kk, ll) * t_F_dX(kk, ll, mm, nn);
1128 for (
int cc = 0; cc != nbCols / 3; ++cc) {
1130 t_m(
i,
k) += t_row_stress(
i) * (t_invH(
j,
k) * t_col_diff_base(
j));
1131 t_m(
i,
k) += t_row_dX_stress(
i,
k,
l) * t_col_diff_base(
l);
1132 t_m(
i,
k) += t_row_stress_dX(
i,
k,
l) * t_col_diff_base(
l);
1145 for (; rr != row_nb_base_fun; ++rr)
1160 HookeElement::OpAleLhs_dX_dX<S>::OpAleLhs_dX_dX(
1161 const std::string row_field,
const std::string col_field,
1162 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1163 :
OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
1174 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2), &
m(
r + 1,
c + 0),
1175 &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2), &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1),
1187 double vol = getVolume();
1190 auto t_w = getFTensor0IntegrationWeight();
1194 const int row_nb_base_fun = row_data.
getN().size2();
1200 auto t_cauchy_stress =
1201 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1202 auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
1203 auto t_eshelby_stress =
1204 getFTensor2FromMat<3, 3>(*dataAtPts->eshelbyStressMat);
1205 auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1206 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1207 auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
1208 auto &det_H = *dataAtPts->detHVec;
1211 for (
int gg = 0; gg != nbIntegrationPts; ++gg) {
1214 double a = t_w * vol * det_H[gg];
1217 t_F_dX(
i,
j,
k,
l) = -(t_h(
i,
m) * t_invH(
m,
k)) * t_invH(
l,
j);
1220 t_D_strain_dX(
i,
j,
m,
n) = 0.;
1221 for (
int ii = 0; ii != 3; ++ii)
1222 for (
int jj = 0; jj != 3; ++jj)
1223 for (
int ll = 0; ll != 3; ++ll)
1224 for (
int kk = 0; kk != 3; ++kk) {
1225 auto &v = t_D_strain_dX(ii, jj, kk, ll);
1226 for (
int mm = 0; mm != 3; ++mm)
1227 for (
int nn = 0; nn != 3; ++nn)
1228 v += t_D(ii, jj, mm, nn) * t_F_dX(mm, nn, kk, ll);
1232 t_eshelby_stress_dX(
i,
j,
m,
n) = t_F(
k,
i) * t_D_strain_dX(
k,
j,
m,
n);
1234 for (
int ii = 0; ii != 3; ++ii)
1235 for (
int jj = 0; jj != 3; ++jj)
1236 for (
int mm = 0; mm != 3; ++mm)
1237 for (
int nn = 0; nn != 3; ++nn) {
1238 auto &v = t_eshelby_stress_dX(ii, jj, mm, nn);
1239 for (
int kk = 0; kk != 3; ++kk)
1240 v += t_F_dX(kk, ii, mm, nn) * t_cauchy_stress(kk, jj);
1243 t_eshelby_stress_dX(
i,
j,
k,
l) *= -1;
1246 t_energy_dX(
k,
l) = t_F_dX(
i,
j,
k,
l) * t_cauchy_stress(
i,
j);
1247 t_energy_dX(
k,
l) +=
1248 (t_strain(
m,
n) * t_D(
m,
n,
i,
j)) * t_F_dX(
i,
j,
k,
l);
1249 t_energy_dX(
k,
l) /= 2.;
1251 for (
int kk = 0; kk != 3; ++kk)
1252 for (
int ll = 0; ll != 3; ++ll) {
1253 auto v = t_energy_dX(kk, ll);
1254 for (
int ii = 0; ii != 3; ++ii)
1255 t_eshelby_stress_dX(ii, ii, kk, ll) += v;
1260 for (; rr != nbRows / 3; ++rr) {
1263 auto t_m = get_tensor2(
K, 3 * rr, 0);
1266 t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
1269 t_row_stress(
i) = a * t_row_diff_base_pulled(
j) * t_eshelby_stress(
i,
j);
1272 t_row_diff_base_pulled_dX(
j,
k,
l) =
1273 -(t_row_diff_base(
i) * t_invH(
i,
k)) * t_invH(
l,
j);
1276 t_row_dX_stress(
i,
k,
l) =
1277 a * (t_row_diff_base_pulled_dX(
j,
k,
l) * t_eshelby_stress(
i,
j));
1280 t_row_stress_dX(
i,
m,
n) =
1281 a * t_row_diff_base_pulled(
j) * t_eshelby_stress_dX(
i,
j,
m,
n);
1287 for (
int cc = 0; cc != nbCols / 3; ++cc) {
1289 t_m(
i,
k) += t_row_stress(
i) * (t_invH(
j,
k) * t_col_diff_base(
j));
1290 t_m(
i,
k) += t_row_dX_stress(
i,
k,
l) * t_col_diff_base(
l);
1291 t_m(
i,
k) += t_row_stress_dX(
i,
k,
l) * t_col_diff_base(
l);
1304 for (; rr != row_nb_base_fun; ++rr)
1322 HookeElement::OpAleLhsPre_dX_dx<S>::OpAleLhsPre_dX_dx(
1323 const std::string row_field,
const std::string col_field,
1324 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1326 dataAtPts(data_at_pts) {
1332 EntityType row_type,
1336 const int nb_integration_pts = row_data.
getN().size1();
1338 auto get_eshelby_stress_dx = [
this, nb_integration_pts]() {
1340 t_eshelby_stress_dx;
1341 dataAtPts->eshelbyStress_dx->resize(81, nb_integration_pts,
false);
1343 for (
int ii = 0; ii != 3; ++ii)
1344 for (
int jj = 0; jj != 3; ++jj)
1345 for (
int kk = 0; kk != 3; ++kk)
1346 for (
int ll = 0; ll != 3; ++ll)
1347 t_eshelby_stress_dx.ptr(ii, jj, kk, ll) =
1348 &(*dataAtPts->eshelbyStress_dx)(mm++, 0);
1349 return t_eshelby_stress_dx;
1352 auto t_eshelby_stress_dx = get_eshelby_stress_dx();
1365 auto t_cauchy_stress =
1366 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1367 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1368 auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
1370 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1372 t_eshelby_stress_dx(
i,
j,
m,
n) =
1373 (t_F(
k,
i) * t_D(
k,
j,
m,
l)) * t_invH(
n,
l);
1374 for (
int ii = 0; ii != 3; ++ii)
1375 for (
int jj = 0; jj != 3; ++jj)
1376 for (
int mm = 0; mm != 3; ++mm)
1377 for (
int nn = 0; nn != 3; ++nn) {
1378 auto &v = t_eshelby_stress_dx(ii, jj, mm, nn);
1379 v += t_invH(nn, ii) * t_cauchy_stress(mm, jj);
1381 t_eshelby_stress_dx(
i,
j,
k,
l) *= -1;
1384 t_energy_dx(
m,
n) = t_invH(
n,
j) * t_cauchy_stress(
m,
j);
1386 for (
int mm = 0; mm != 3; ++mm)
1387 for (
int nn = 0; nn != 3; ++nn) {
1388 auto v = t_energy_dx(mm, nn);
1389 for (
int ii = 0; ii != 3; ++ii)
1390 t_eshelby_stress_dx(ii, ii, mm, nn) += v;
1396 ++t_eshelby_stress_dx;
1403 template <
class ELEMENT>
1404 HookeElement::OpPostProcHookeElement<ELEMENT>::OpPostProcHookeElement(
1405 const string row_field, boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1407 std::vector<EntityHandle> &map_gauss_pts,
bool is_ale,
bool is_field_disp)
1409 dataAtPts(data_at_pts), blockSetsPtr(block_sets_ptr),
1410 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts), isALE(is_ale),
1411 isFieldDisp(is_field_disp) {}
1413 template <
class ELEMENT>
1418 if (
type != MBVERTEX) {
1422 auto tensor_to_tensor = [](
const auto &t1,
auto &t2) {
1423 t2(0, 0) = t1(0, 0);
1424 t2(1, 1) = t1(1, 1);
1425 t2(2, 2) = t1(2, 2);
1426 t2(0, 1) = t2(1, 0) = t1(1, 0);
1427 t2(0, 2) = t2(2, 0) = t1(2, 0);
1428 t2(1, 2) = t2(2, 1) = t1(2, 1);
1431 std::array<double,9> def_val;
1434 auto make_tag = [&](
auto name,
auto size) {
1436 CHKERR postProcMesh.tag_get_handle(name, size, MB_TYPE_DOUBLE, th,
1437 MB_TAG_CREAT | MB_TAG_SPARSE,
1442 auto th_stress = make_tag(
"STRESS", 9);
1443 auto th_psi = make_tag(
"ENERGY", 1);
1445 const int nb_integration_pts = mapGaussPts.size();
1452 auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1453 auto t_H = getFTensor2FromMat<3, 3>(*dataAtPts->HMat);
1455 dataAtPts->stiffnessMat->resize(36, 1,
false);
1458 for (
auto &
m : (blockSetsPtr)) {
1459 const double young =
m.second.E;
1460 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;
1483 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1492 t_small_strain_symm(
i,
j) = (t_h(
i,
j) || t_h(
j,
i)) / 2.;
1496 t_F(
i,
j) = t_h(
i,
k) * t_invH(
k,
j);
1497 t_small_strain_symm(
i,
j) = (t_F(
i,
j) || t_F(
j,
i)) / 2.;
1501 t_small_strain_symm(0, 0) -= 1;
1502 t_small_strain_symm(1, 1) -= 1;
1503 t_small_strain_symm(2, 2) -= 1;
1506 t_stress_symm(
i,
j) = t_D(
i,
j,
k,
l) * t_small_strain_symm(
k,
l);
1507 tensor_to_tensor(t_stress_symm, t_stress);
1509 const double psi = 0.5 * t_stress_symm(
i,
j) * t_small_strain_symm(
i,
j);
1511 CHKERR postProcMesh.tag_set_data(th_psi, &mapGaussPts[gg], 1, &psi);
1512 CHKERR postProcMesh.tag_set_data(th_stress, &mapGaussPts[gg], 1,
1522 HookeElement::OpAnalyticalInternalStain_dx<S>::OpAnalyticalInternalStain_dx(
1523 const std::string row_field,
1524 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
1526 :
OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1527 strainFun(strain_fun) {}
1531 HookeElement::OpAnalyticalInternalStain_dx<S>::iNtegrate(
EntData &row_data) {
1540 &v(
r + 0), &v(
r + 1), &v(
r + 2));
1543 const int nb_integration_pts = getGaussPts().size2();
1545 auto get_coords = [&]() {
1546 if (getHoCoordsAtGaussPts().size1() == nb_integration_pts)
1547 return getFTensor1HoCoordsAtGaussPts();
1549 return getFTensor1CoordsAtGaussPts();
1551 auto t_coords = get_coords();
1554 double vol = getVolume();
1555 auto t_w = getFTensor0IntegrationWeight();
1557 nF.resize(nbRows,
false);
1567 const int row_nb_base_fun = row_data.
getN().size2();
1569 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
1571 auto t_fun_strain = strainFun(t_coords);
1573 t_stress(
i,
j) = t_D(
i,
j,
k,
l) * t_fun_strain(
k,
l);
1576 double a = t_w * vol;
1578 if (getHoGaussPtsDetJac().size()) {
1580 a *= getHoGaussPtsDetJac()[gg];
1583 auto t_nf = get_tensor1(nF, 0);
1586 for (; rr != nbRows / 3; ++rr) {
1587 t_nf(
i) += a * t_row_diff_base(
j) * t_stress(
i,
j);
1592 for (; rr != row_nb_base_fun; ++rr)
1604 HookeElement::OpAnalyticalInternalAleStain_dX<S>::
1605 OpAnalyticalInternalAleStain_dX(
1606 const std::string row_field,
1607 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
1609 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr)
1610 :
OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1611 strainFun(strain_fun), matPosAtPtsPtr(mat_pos_at_pts_ptr) {}
1615 HookeElement::OpAnalyticalInternalAleStain_dX<S>::iNtegrate(
EntData &row_data) {
1624 &v(
r + 0), &v(
r + 1), &v(
r + 2));
1627 const int nb_integration_pts = getGaussPts().size2();
1629 auto get_coords = [&]() {
return getFTensor1FromMat<3>(*matPosAtPtsPtr); };
1630 auto t_coords = get_coords();
1633 double vol = getVolume();
1634 auto t_w = getFTensor0IntegrationWeight();
1636 nF.resize(nbRows,
false);
1643 auto t_F = getFTensor2FromMat<3, 3>(*(dataAtPts->FMat));
1644 auto &det_H = *dataAtPts->detHVec;
1645 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1649 const int row_nb_base_fun = row_data.
getN().size2();
1651 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
1653 auto t_fun_strain = strainFun(t_coords);
1655 t_stress(
i,
j) = t_D(
i,
j,
k,
l) * t_fun_strain(
k,
l);
1657 t_eshelby_stress(
i,
j) = -t_F(
k,
i) * t_stress(
k,
j);
1660 double a = t_w * vol * det_H[gg];
1662 auto t_nf = get_tensor1(nF, 0);
1665 for (; rr != nbRows / 3; ++rr) {
1667 t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
1668 t_nf(
i) += a * t_row_diff_base_pulled(
j) * t_eshelby_stress(
i,
j);
1673 for (; rr != row_nb_base_fun; ++rr)
1687 HookeElement::OpAnalyticalInternalAleStain_dx<S>::
1688 OpAnalyticalInternalAleStain_dx(
1689 const std::string row_field,
1690 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
1692 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr)
1693 :
OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1694 strainFun(strain_fun), matPosAtPtsPtr(mat_pos_at_pts_ptr) {}
1698 HookeElement::OpAnalyticalInternalAleStain_dx<S>::iNtegrate(
EntData &row_data) {
1707 &v(
r + 0), &v(
r + 1), &v(
r + 2));
1710 const int nb_integration_pts = getGaussPts().size2();
1712 auto get_coords = [&]() {
return getFTensor1FromMat<3>(*matPosAtPtsPtr); };
1713 auto t_coords = get_coords();
1716 double vol = getVolume();
1717 auto t_w = getFTensor0IntegrationWeight();
1719 nF.resize(nbRows,
false);
1726 auto &det_H = *dataAtPts->detHVec;
1727 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1731 const int row_nb_base_fun = row_data.
getN().size2();
1733 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
1735 auto t_fun_strain = strainFun(t_coords);
1737 t_stress(
i,
j) = t_D(
i,
j,
k,
l) * t_fun_strain(
k,
l);
1740 double a = t_w * vol * det_H[gg];
1742 auto t_nf = get_tensor1(nF, 0);
1745 for (; rr != nbRows / 3; ++rr) {
1747 t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
1748 t_nf(
i) += a * t_row_diff_base_pulled(
j) * t_stress(
i,
j);
1753 for (; rr != row_nb_base_fun; ++rr)
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, Vec *v_energy_ptr)
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)
DataForcesAndSourcesCore::EntData EntData
ForcesAndSourcesCore::UserDataOperator UserDataOperator
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)
#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 ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
struct ConvectiveMassElement BlockData
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
FTensor::Index< 'k', 3 > k
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
ublas::vector< int, IntAllocator > VectorInt
ublas::vector< double, DoubleAllocator > VectorDouble
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
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.
DeprecatedCoreInterface Interface
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
const double D
diffusivity
const double r
rate factor
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
Data on single entity (This is passed as argument to DataOperator::doWork)
const VectorInt & getIndices() const
Get global indices of dofs on entity.
FTensor::Tensor1< double *, 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....
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Deprecated interface functions.
Data operator to do calculations at integration points.
default operator for TET element
Volume finite element with switches.
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)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
boost::function< FTensor::Tensor2_symmetric< double, 3 > FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &t_coords) > StrainFunctions
StrainFunctions strainFun
boost::shared_ptr< MatrixDouble > matPosAtPtsPtr
StrainFunctions strainFun
boost::function< FTensor::Tensor2_symmetric< double, 3 > FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &t_coords) > StrainFunctions
boost::shared_ptr< MatrixDouble > matPosAtPtsPtr
boost::function< FTensor::Tensor2_symmetric< double, 3 > FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, 3 > &t_coords) > StrainFunctions
StrainFunctions strainFun
int nbCols
number if dof on column
bool isDiag
true if this block is on diagonal
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
int nbIntegrationPts
number of integration points
int nbRows
number of dofs on rows
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< map< int, BlockData > > blockSetsPtr
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, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
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
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
const double rhoN
exponent n in E(p) = E * (p / p_0)^n
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
std::vector< EntityHandle > & mapGaussPts
map< int, BlockData > & blockSetsPtr
moab::Interface & postProcMesh
boost::shared_ptr< DataAtIntegrationPts > dataAtPts