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>>
213 :
public VolumeElementForcesAndSourcesCore::UserDataOperator {
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) {}
231 PetscErrorCode
doWork(
int row_side,
int col_side, EntityType row_type,
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);
356 MoFEMErrorCode doWork(
int row_side,
int col_side, EntityType row_type,
357 EntityType col_type,
EntData &row_data,
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);
436 OpAleLhs_dx_dx(
const std::string row_field,
const std::string col_field,
437 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
451 OpAleLhs_dx_dX(
const std::string row_field,
const std::string col_field,
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);
523 OpAleLhs_dX_dX(
const std::string row_field,
const std::string col_field,
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,
673 const EntityType type = MBTET,
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) {
703 const int nb_integration_pts = getGaussPts().size2();
704 dataAtPts->smallStrainMat->resize(6, nb_integration_pts,
false);
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) {}
737 auto get_tensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
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();
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) {
837 int row_side, EntityType row_type,
EntData &row_data) {
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) {
898 const int nb_integration_pts = getGaussPts().size2();
900 dataAtPts->cauchyStressMat->resize(6, nb_integration_pts,
false);
901 auto t_cauchy_stress =
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) {}
935 auto get_tensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
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) {}
1025 auto get_tensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
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 =
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) {}
1154 auto get_tensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
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 =
1185 auto t_eshelby_stress =
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) {
1314 EntityType row_type,
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 =
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>
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();
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;
1459 for (
auto &
m : (blockSetsPtr)) {
1460 if (
m.second.tEts.find(ent_3d) !=
m.second.tEts.end()) {
1461 const double young =
m.second.E;
1462 const double poisson =
m.second.PoissonRatio;
1463 const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
1464 block_id =
m.second.iD;
1466 t_D(
i,
j,
k,
l) = 0.;
1467 t_D(0, 0, 0, 0) = t_D(1, 1, 1, 1) = t_D(2, 2, 2, 2) = 1 - poisson;
1468 t_D(0, 1, 0, 1) = t_D(0, 2, 0, 2) = t_D(1, 2, 1, 2) =
1469 0.5 * (1 - 2 * poisson);
1470 t_D(0, 0, 1, 1) = t_D(1, 1, 0, 0) = t_D(0, 0, 2, 2) = t_D(2, 2, 0, 0) =
1471 t_D(1, 1, 2, 2) = t_D(2, 2, 1, 1) = poisson;
1472 t_D(
i,
j,
k,
l) *= coefficient;
1480 "Element not found in any of material blocksets");
1482 int def_val_int = 0;
1484 CHKERR postProcMesh.tag_get_handle(
"MAT_ELASTIC", 1, MB_TYPE_INTEGER, tag_mat,
1485 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val_int);
1494 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1503 t_small_strain_symm(
i,
j) = (t_h(
i,
j) || t_h(
j,
i)) / 2.;
1507 t_F(
i,
j) = t_h(
i,
k) * t_invH(
k,
j);
1508 t_small_strain_symm(
i,
j) = (t_F(
i,
j) || t_F(
j,
i)) / 2.;
1512 t_small_strain_symm(0, 0) -= 1;
1513 t_small_strain_symm(1, 1) -= 1;
1514 t_small_strain_symm(2, 2) -= 1;
1517 t_stress_symm(
i,
j) = t_D(
i,
j,
k,
l) * t_small_strain_symm(
k,
l);
1518 tensor_to_tensor(t_stress_symm, t_stress);
1520 const double psi = 0.5 * t_stress_symm(
i,
j) * t_small_strain_symm(
i,
j);
1522 CHKERR postProcMesh.tag_set_data(th_psi, &mapGaussPts[gg], 1, &psi);
1523 CHKERR postProcMesh.tag_set_data(th_stress, &mapGaussPts[gg], 1,
1525 CHKERR postProcMesh.tag_set_data(tag_mat, &mapGaussPts[gg], 1,
1535HookeElement::OpAnalyticalInternalStrain_dx<S>::OpAnalyticalInternalStrain_dx(
1536 const std::string row_field,
1537 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1539 :
OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1540 strainFun(strain_fun) {}
1544HookeElement::OpAnalyticalInternalStrain_dx<S>::iNtegrate(
EntData &row_data) {
1553 &
v(r + 0), &
v(r + 1), &
v(r + 2));
1556 const int nb_integration_pts = getGaussPts().size2();
1557 auto t_coords = getFTensor1CoordsAtGaussPts();
1560 double vol = getVolume();
1561 auto t_w = getFTensor0IntegrationWeight();
1563 nF.resize(nbRows,
false);
1573 const int row_nb_base_fun = row_data.
getN().size2();
1575 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
1577 auto t_fun_strain = strainFun(t_coords);
1579 t_stress(
i,
j) = -t_D(
i,
j,
k,
l) * t_fun_strain(
k,
l);
1582 double a = t_w * vol;
1584 auto t_nf = get_tensor1(nF, 0);
1587 for (; rr != nbRows / 3; ++rr) {
1588 t_nf(
i) +=
a * t_row_diff_base(
j) * t_stress(
i,
j);
1593 for (; rr != row_nb_base_fun; ++rr)
1605HookeElement::OpAnalyticalInternalAleStrain_dX<S>::
1606 OpAnalyticalInternalAleStrain_dX(
1607 const std::string row_field,
1608 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
1610 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr)
1611 :
OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1612 strainFun(strain_fun), matPosAtPtsPtr(mat_pos_at_pts_ptr) {}
1625 &
v(r + 0), &
v(r + 1), &
v(r + 2));
1628 const int nb_integration_pts = getGaussPts().size2();
1631 auto t_coords = get_coords();
1634 double vol = getVolume();
1635 auto t_w = getFTensor0IntegrationWeight();
1637 nF.resize(nbRows,
false);
1645 auto &det_H = *dataAtPts->detHVec;
1650 const int row_nb_base_fun = row_data.
getN().size2();
1652 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
1654 auto t_fun_strain = strainFun(t_coords);
1656 t_stress(
i,
j) = -t_D(
i,
j,
k,
l) * t_fun_strain(
k,
l);
1658 t_eshelby_stress(
i,
j) = -t_F(
k,
i) * t_stress(
k,
j);
1661 double a = t_w * vol * det_H[gg];
1663 auto t_nf = get_tensor1(nF, 0);
1666 for (; rr != nbRows / 3; ++rr) {
1668 t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
1669 t_nf(
i) +=
a * t_row_diff_base_pulled(
j) * t_eshelby_stress(
i,
j);
1674 for (; rr != row_nb_base_fun; ++rr)
1688HookeElement::OpAnalyticalInternalAleStrain_dx<S>::
1689 OpAnalyticalInternalAleStrain_dx(
1690 const std::string row_field,
1691 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
1693 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr)
1694 :
OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1695 strainFun(strain_fun), matPosAtPtsPtr(mat_pos_at_pts_ptr) {}
1708 &
v(r + 0), &
v(r + 1), &
v(r + 2));
1711 const int nb_integration_pts = getGaussPts().size2();
1714 auto t_coords = get_coords();
1717 double vol = getVolume();
1718 auto t_w = getFTensor0IntegrationWeight();
1720 nF.resize(nbRows,
false);
1727 auto &det_H = *dataAtPts->detHVec;
1732 const int row_nb_base_fun = row_data.
getN().size2();
1734 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
1736 auto t_fun_strain = strainFun(t_coords);
1738 t_stress(
i,
j) = -t_D(
i,
j,
k,
l) * t_fun_strain(
k,
l);
1741 double a = t_w * vol * det_H[gg];
1743 auto t_nf = get_tensor1(nF, 0);
1746 for (; rr != nbRows / 3; ++rr) {
1748 t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
1749 t_nf(
i) +=
a * t_row_diff_base_pulled(
j) * t_stress(
i,
j);
1754 for (; rr != row_nb_base_fun; ++rr)
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 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 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< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
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.
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
static auto getFTensor2SymmetricFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 2 (matrix) form data matrix.
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.
FTensor::Index< 'm', 3 > m
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
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.
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
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
OpAleLhs_dX_dx(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
boost::shared_ptr< MatrixDouble > matPosAtPtsPtr
boost::function< FTensor::Tensor2_symmetric< double, 3 >(FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &t_coords) > StrainFunction
boost::shared_ptr< MatrixDouble > matPosAtPtsPtr
boost::function< FTensor::Tensor2_symmetric< double, 3 >(FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &t_coords) > StrainFunction
boost::function< FTensor::Tensor2_symmetric< double, 3 >(FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, 3 > &t_coords) > StrainFunction
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
SmartPetscObj< Vec > ghostVec
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
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, 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
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