Operators and data structures for linear elastic analysis.
In other words, spatial deformation is small but topological changes large.
#ifndef __HOOKE_ELEMENT_HPP
#define __HOOKE_ELEMENT_HPP
#ifndef __BASICFINITEELEMENTS_HPP__
#endif
#ifndef __NONLINEAR_ELASTIC_HPP
};
};
#endif
#ifndef __CONVECTIVE_MASS_ELEMENT_HPP
};
}
#endif
using EntData = EntitiesFieldData::EntData;
VolumeElementForcesAndSourcesCore::UserDataOperator;
boost::shared_ptr<MatrixDouble> smallStrainMat;
boost::shared_ptr<MatrixDouble> hMat;
boost::shared_ptr<MatrixDouble> FMat;
boost::shared_ptr<MatrixDouble> HMat;
boost::shared_ptr<VectorDouble> detHVec;
boost::shared_ptr<MatrixDouble> invHMat;
boost::shared_ptr<MatrixDouble> cauchyStressMat;
boost::shared_ptr<MatrixDouble> stiffnessMat;
boost::shared_ptr<VectorDouble> energyVec;
boost::shared_ptr<MatrixDouble> eshelbyStressMat;
boost::shared_ptr<MatrixDouble> eshelbyStress_dx;
smallStrainMat = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
detHVec = boost::shared_ptr<VectorDouble>(
new VectorDouble());
invHMat = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
cauchyStressMat = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
stiffnessMat = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
energyVec = boost::shared_ptr<VectorDouble>(
new VectorDouble());
eshelbyStressMat = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
eshelbyStress_dx = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
}
Range forcesOnlyOnEntitiesRow;
Range forcesOnlyOnEntitiesCol;
};
template <bool D = false>
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
private:
boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
};
const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
private:
boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
};
#define MAT_TO_DDG(SM) \
&(*SM)(0, 0), &(*SM)(1, 0), &(*SM)(2, 0), &(*SM)(3, 0), &(*SM)(4, 0), \
&(*SM)(5, 0), &(*SM)(6, 0), &(*SM)(7, 0), &(*SM)(8, 0), &(*SM)(9, 0), \
&(*SM)(10, 0), &(*SM)(11, 0), &(*SM)(12, 0), &(*SM)(13, 0), \
&(*SM)(14, 0), &(*SM)(15, 0), &(*SM)(16, 0), &(*SM)(17, 0), \
&(*SM)(18, 0), &(*SM)(19, 0), &(*SM)(20, 0), &(*SM)(21, 0), \
&(*SM)(22, 0), &(*SM)(23, 0), &(*SM)(24, 0), &(*SM)(25, 0), \
&(*SM)(26, 0), &(*SM)(27, 0), &(*SM)(28, 0), &(*SM)(29, 0), \
&(*SM)(30, 0), &(*SM)(31, 0), &(*SM)(32, 0), &(*SM)(33, 0), \
&(*SM)(34, 0), &(*SM)(35, 0)
boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
protected:
boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
};
boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
SmartPetscObj<Vec> ghost_vec = SmartPetscObj<Vec>());
protected:
boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
SmartPetscObj<Vec> ghostVec;
};
const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
protected:
boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
};
template <int S = 0>
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);
protected:
boost::shared_ptr<map<int, BlockData>>
blockSetsPtr;
boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
};
: public VolumeElementForcesAndSourcesCore::UserDataOperator {
boost::shared_ptr<DataAtIntegrationPts> commonData;
const std::string col_field,
BlockData &data,
boost::shared_ptr<DataAtIntegrationPts> &common_data,
bool symm = true)
row_field, col_field, OPROWCOL, symm),
commonData(common_data), dAta(data), massData(mass_data) {}
PetscErrorCode doWork(
int row_side,
int col_side,
EntityType row_type,
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
&
m(3 * r + 0, 3 *
c + 0), &
m(3 * r + 0, 3 *
c + 1),
&
m(3 * r + 0, 3 *
c + 2), &
m(3 * r + 1, 3 *
c + 0),
&
m(3 * r + 1, 3 *
c + 1), &
m(3 * r + 1, 3 *
c + 2),
&
m(3 * r + 2, 3 *
c + 0), &
m(3 * r + 2, 3 *
c + 1),
&
m(3 * r + 2, 3 *
c + 2));
};
const int row_nb_dofs = row_data.getIndices().size();
if (!row_nb_dofs)
const int col_nb_dofs = col_data.getIndices().size();
if (!col_nb_dofs)
if (dAta.
tEts.find(getFEEntityHandle()) == dAta.
tEts.end()) {
}
if (massData.
tEts.find(getFEEntityHandle()) == massData.
tEts.end()) {
}
const bool diagonal_block =
(row_type == col_type) && (row_side == col_side);
locK.resize(row_nb_dofs, col_nb_dofs, false);
locK.clear();
const int row_nb_gauss_pts = row_data.getN().size1();
const int row_nb_base_functions = row_data.getN().size2();
double density = massData.
rho0;
auto t_w = getFTensor0IntegrationWeight();
for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
auto t_row_base_func = row_data.getFTensor0N(gg, 0);
double w = getVolume() * t_w;
for (int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
auto t_col_base_func = col_data.getFTensor0N(gg, 0);
for (int col_bb = 0; col_bb != col_nb_dofs / 3; col_bb++) {
auto t_assemble = get_tensor2(locK, row_bb, col_bb);
t_assemble(
i,
j) += density * t_row_base_func * t_col_base_func *
w;
++t_col_base_func;
}
++t_row_base_func;
}
++t_w;
}
ADD_VALUES);
if (row_type != col_type || row_side != col_side) {
translocK.resize(col_nb_dofs, row_nb_dofs, false);
noalias(translocK) = trans(locK);
ADD_VALUES);
}
}
};
protected:
boost::shared_ptr<map<int, BlockData>>
blockSetsPtr;
boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
const double rhoN;
const double rHo0;
public:
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);
};
OpAssemble(
const std::string row_field,
const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
const char type, bool symm = false);
protected:
boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
int nbRows;
int nbCols;
int nbIntegrationPts;
bool isDiag;
};
OpRhs_dx(
const std::string row_field,
const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
protected:
};
OpLhs_dx_dx(
const std::string row_field,
const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
protected:
};
OpAleRhs_dx(
const std::string row_field,
const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
protected:
};
OpAleLhs_dx_dx(
const std::string row_field,
const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
protected:
};
OpAleLhs_dx_dX(
const std::string row_field,
const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
protected:
};
boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
boost::shared_ptr<MatrixDouble> rhoGradAtGaussPtsPtr;
const double rhoN;
const double rHo0;
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);
protected:
};
boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
boost::shared_ptr<MatrixDouble> rhoGradAtGaussPtsPtr;
const double rhoN;
const double rHo0;
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);
protected:
};
OpAleRhs_dX(
const std::string row_field,
const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
protected:
};
OpAleLhs_dX_dX(
const std::string row_field,
const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
protected:
};
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
private:
boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
};
OpAleLhs_dX_dx(
const std::string row_field,
const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
:
OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false) {}
protected:
};
typedef boost::function<
)
>
StrainFunction;
const std::string row_field,
boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
StrainFunction strain_fun);
protected:
StrainFunction strainFun;
};
typedef boost::function<
)
>
StrainFunction;
const std::string row_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
StrainFunction strain_fun,
boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr);
protected:
StrainFunction strainFun;
boost::shared_ptr<MatrixDouble> matPosAtPtsPtr;
};
typedef boost::function<
)
>
StrainFunction;
const std::string row_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
StrainFunction strain_fun,
boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr);
protected:
StrainFunction strainFun;
boost::shared_ptr<MatrixDouble> matPosAtPtsPtr;
};
template <class ELEMENT>
boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
map<int, BlockData>
&blockSetsPtr;
moab::Interface &postProcMesh;
std::vector<EntityHandle> &mapGaussPts;
bool isALE;
bool isFieldDisp;
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);
EntitiesFieldData::EntData &data);
};
boost::shared_ptr<map<int, BlockData>> &block_sets_ptr);
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);
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,
boost::shared_ptr<DataAtIntegrationPts> data_at_pts = nullptr);
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);
private:
};
template <bool D>
HookeElement::OpCalculateStrain<D>::OpCalculateStrain(
const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
dataAtPts(data_at_pts) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
template <bool D>
MoFEMErrorCode HookeElement::OpCalculateStrain<D>::doWork(
int row_side,
const int nb_integration_pts = getGaussPts().size2();
dataAtPts->smallStrainMat->resize(6, nb_integration_pts, false);
auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
auto t_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hMat));
for (int gg = 0; gg != nb_integration_pts; ++gg) {
t_strain(
i,
j) = (t_h(
i,
j) || t_h(
j,
i)) / 2.;
t_strain(0, 0) -= 1;
t_strain(1, 1) -= 1;
t_strain(2, 2) -= 1;
}
++t_strain;
++t_h;
}
}
template <int S>
HookeElement::OpAleLhs_dx_dx<S>::OpAleLhs_dx_dx(
const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
:
OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
template <int S>
&
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 1,
c + 0),
&
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 2,
c + 0), &
m(r + 2,
c + 1),
};
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
const int row_nb_base_fun = row_data.
getN().size2();
auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
auto &det_H = *dataAtPts->detHVec;
for (int gg = 0; gg != nbIntegrationPts; ++gg) {
double a = t_w * vol * det_H[gg];
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
auto t_m = get_tensor2(
K, 3 * rr, 0);
t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
t_rowD(
l,
j,
k) = t_D(
i,
j,
k,
l) * (
a * t_row_diff_base_pulled(
i));
for (int cc = 0; cc != nbCols / 3; ++cc) {
t_col_diff_base_pulled(
j) = t_col_diff_base(
i) * t_invH(
i,
j);
t_m(
i,
j) += t_rowD(
i,
j,
k) * t_col_diff_base_pulled(
k);
++t_col_diff_base;
++t_m;
}
++t_row_diff_base;
}
for (; rr != row_nb_base_fun; ++rr)
++t_row_diff_base;
++t_w;
++t_D;
++t_invH;
}
}
template <int S>
HookeElement::OpCalculateHomogeneousStiffness<S>::
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)
blockSetsPtr(block_sets_ptr), dataAtPts(data_at_pts) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
template <int S>
MoFEMErrorCode HookeElement::OpCalculateHomogeneousStiffness<S>::doWork(
for (
auto &
m : (*blockSetsPtr)) {
if (
m.second.tEts.find(getFEEntityHandle()) !=
m.second.tEts.end()) {
dataAtPts->stiffnessMat->resize(36, 1, false);
const double young =
m.second.E;
const double poisson =
m.second.PoissonRatio;
const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
t_D(0, 0, 0, 0) = 1 - poisson;
t_D(1, 1, 1, 1) = 1 - poisson;
t_D(2, 2, 2, 2) = 1 - poisson;
t_D(0, 1, 0, 1) = 0.5 * (1 - 2 * poisson);
t_D(0, 2, 0, 2) = 0.5 * (1 - 2 * poisson);
t_D(1, 2, 1, 2) = 0.5 * (1 - 2 * poisson);
t_D(0, 0, 1, 1) = poisson;
t_D(1, 1, 0, 0) = poisson;
t_D(0, 0, 2, 2) = poisson;
t_D(2, 2, 0, 0) = poisson;
t_D(1, 1, 2, 2) = poisson;
t_D(2, 2, 1, 1) = poisson;
t_D(
i,
j,
k,
l) *= coefficient;
break;
}
}
}
template <int S>
HookeElement::OpCalculateStress<S>::OpCalculateStress(
const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
dataAtPts(data_at_pts) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
template <int S>
MoFEMErrorCode HookeElement::OpCalculateStress<S>::doWork(
int row_side,
const int nb_integration_pts = getGaussPts().size2();
auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
dataAtPts->cauchyStressMat->resize(6, nb_integration_pts, false);
auto t_cauchy_stress =
getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
for (int gg = 0; gg != nb_integration_pts; ++gg) {
t_cauchy_stress(
i,
j) = t_D(
i,
j,
k,
l) * t_strain(
k,
l);
++t_strain;
++t_cauchy_stress;
++t_D;
}
}
template <int S>
HookeElement::OpLhs_dx_dx<S>::OpLhs_dx_dx(
const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
:
OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
template <int S>
&
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 1,
c + 0),
&
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 2,
c + 0), &
m(r + 2,
c + 1),
};
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
const int row_nb_base_fun = row_data.
getN().size2();
for (int gg = 0; gg != nbIntegrationPts; ++gg) {
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
auto t_m = get_tensor2(
K, 3 * rr, 0);
t_rowD(
l,
j,
k) = t_D(
i,
j,
k,
l) * (
a * t_row_diff_base(
i));
for (int cc = 0; cc != nbCols / 3; ++cc) {
t_m(
i,
j) += t_rowD(
i,
j,
k) * t_col_diff_base(
k);
++t_col_diff_base;
++t_m;
}
++t_row_diff_base;
}
for (; rr != row_nb_base_fun; ++rr)
++t_row_diff_base;
++t_w;
++t_D;
}
}
template <int S>
HookeElement::OpAleLhs_dx_dX<S>::OpAleLhs_dx_dX(
const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
:
OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false) {}
template <int S>
&
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 1,
c + 0),
&
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 2,
c + 0), &
m(r + 2,
c + 1),
};
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
const int row_nb_base_fun = row_data.
getN().size2();
auto t_cauchy_stress =
getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
auto &det_H = *dataAtPts->detHVec;
for (int gg = 0; gg != nbIntegrationPts; ++gg) {
double a = t_w * vol * det_H[gg];
t_F_dX(
i,
j,
k,
l) = -(t_h(
i,
m) * t_invH(
m,
k)) * t_invH(
l,
j);
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
auto t_m = get_tensor2(
K, 3 * rr, 0);
t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
t_row_stress(
i) =
a * t_row_diff_base_pulled(
j) * t_cauchy_stress(
i,
j);
t_row_diff_base_pulled_dX(
j,
k,
l) =
-(t_invH(
i,
k) * t_row_diff_base(
i)) * t_invH(
l,
j);
t_row_dX_stress(
i,
k,
l) =
a * (t_row_diff_base_pulled_dX(
j,
k,
l) * t_cauchy_stress(
j,
i));
t_row_D(
l,
j,
k) = (
a * t_row_diff_base_pulled(
i)) * t_D(
i,
j,
k,
l);
t_row_stress_dX(
i,
j,
k) = 0;
for (int ii = 0; ii != 3; ++ii)
for (int mm = 0; mm != 3; ++mm)
for (int nn = 0; nn != 3; ++nn) {
auto &
v = t_row_stress_dX(ii, mm, nn);
for (int kk = 0; kk != 3; ++kk)
for (int ll = 0; ll != 3; ++ll)
v += t_row_D(ii, kk, ll) * t_F_dX(kk, ll, mm, nn);
}
for (int cc = 0; cc != nbCols / 3; ++cc) {
t_m(
i,
k) += t_row_stress(
i) * (t_invH(
j,
k) * t_col_diff_base(
j));
t_m(
i,
k) += t_row_dX_stress(
i,
k,
l) * t_col_diff_base(
l);
t_m(
i,
k) += t_row_stress_dX(
i,
k,
l) * t_col_diff_base(
l);
++t_col_diff_base;
++t_m;
}
++t_row_diff_base;
}
for (; rr != row_nb_base_fun; ++rr)
++t_row_diff_base;
++t_w;
++t_D;
++t_cauchy_stress;
++t_invH;
++t_h;
}
}
template <int S>
HookeElement::OpAleLhs_dX_dX<S>::OpAleLhs_dX_dX(
const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
:
OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
template <int S>
&
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 1,
c + 0),
&
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 2,
c + 0), &
m(r + 2,
c + 1),
};
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
const int row_nb_base_fun = row_data.
getN().size2();
auto t_cauchy_stress =
getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
auto t_eshelby_stress =
getFTensor2FromMat<3, 3>(*dataAtPts->eshelbyStressMat);
auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
auto &det_H = *dataAtPts->detHVec;
for (int gg = 0; gg != nbIntegrationPts; ++gg) {
double a = t_w * vol * det_H[gg];
t_F_dX(
i,
j,
k,
l) = -(t_h(
i,
m) * t_invH(
m,
k)) * t_invH(
l,
j);
t_D_strain_dX(
i,
j,
m,
n) = 0.;
for (int ii = 0; ii != 3; ++ii)
for (int jj = 0; jj != 3; ++jj)
for (int ll = 0; ll != 3; ++ll)
for (int kk = 0; kk != 3; ++kk) {
auto &
v = t_D_strain_dX(ii, jj, kk, ll);
for (int mm = 0; mm != 3; ++mm)
for (int nn = 0; nn != 3; ++nn)
v += t_D(ii, jj, mm, nn) * t_F_dX(mm, nn, kk, ll);
}
t_eshelby_stress_dX(
i,
j,
m,
n) = t_F(
k,
i) * t_D_strain_dX(
k,
j,
m,
n);
for (int ii = 0; ii != 3; ++ii)
for (int jj = 0; jj != 3; ++jj)
for (int mm = 0; mm != 3; ++mm)
for (int nn = 0; nn != 3; ++nn) {
auto &
v = t_eshelby_stress_dX(ii, jj, mm, nn);
for (int kk = 0; kk != 3; ++kk)
v += t_F_dX(kk, ii, mm, nn) * t_cauchy_stress(kk, jj);
}
t_eshelby_stress_dX(
i,
j,
k,
l) *= -1;
t_energy_dX(
k,
l) = t_F_dX(
i,
j,
k,
l) * t_cauchy_stress(
i,
j);
(t_strain(
m,
n) * t_D(
m,
n,
i,
j)) * t_F_dX(
i,
j,
k,
l);
for (int kk = 0; kk != 3; ++kk)
for (int ll = 0; ll != 3; ++ll) {
auto v = t_energy_dX(kk, ll);
for (int ii = 0; ii != 3; ++ii)
t_eshelby_stress_dX(ii, ii, kk, ll) +=
v;
}
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
auto t_m = get_tensor2(
K, 3 * rr, 0);
t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
t_row_stress(
i) =
a * t_row_diff_base_pulled(
j) * t_eshelby_stress(
i,
j);
t_row_diff_base_pulled_dX(
j,
k,
l) =
-(t_row_diff_base(
i) * t_invH(
i,
k)) * t_invH(
l,
j);
t_row_dX_stress(
i,
k,
l) =
a * (t_row_diff_base_pulled_dX(
j,
k,
l) * t_eshelby_stress(
i,
j));
t_row_stress_dX(
i,
m,
n) =
a * t_row_diff_base_pulled(
j) * t_eshelby_stress_dX(
i,
j,
m,
n);
for (int cc = 0; cc != nbCols / 3; ++cc) {
t_m(
i,
k) += t_row_stress(
i) * (t_invH(
j,
k) * t_col_diff_base(
j));
t_m(
i,
k) += t_row_dX_stress(
i,
k,
l) * t_col_diff_base(
l);
t_m(
i,
k) += t_row_stress_dX(
i,
k,
l) * t_col_diff_base(
l);
++t_col_diff_base;
++t_m;
}
++t_row_diff_base;
}
for (; rr != row_nb_base_fun; ++rr)
++t_row_diff_base;
++t_w;
++t_D;
++t_cauchy_stress;
++t_strain;
++t_eshelby_stress;
++t_h;
++t_invH;
++t_F;
}
}
template <int S>
HookeElement::OpAleLhsPre_dX_dx<S>::OpAleLhsPre_dX_dx(
const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
dataAtPts(data_at_pts) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
template <int S>
MoFEMErrorCode HookeElement::OpAleLhsPre_dX_dx<S>::doWork(
int row_side,
const int nb_integration_pts = row_data.
getN().size1();
auto get_eshelby_stress_dx = [this, nb_integration_pts]() {
t_eshelby_stress_dx;
dataAtPts->eshelbyStress_dx->resize(81, nb_integration_pts, false);
int mm = 0;
for (int ii = 0; ii != 3; ++ii)
for (int jj = 0; jj != 3; ++jj)
for (int kk = 0; kk != 3; ++kk)
for (int ll = 0; ll != 3; ++ll)
t_eshelby_stress_dx.ptr(ii, jj, kk, ll) =
&(*dataAtPts->eshelbyStress_dx)(mm++, 0);
return t_eshelby_stress_dx;
};
auto t_eshelby_stress_dx = get_eshelby_stress_dx();
auto t_cauchy_stress =
getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
for (int gg = 0; gg != nb_integration_pts; ++gg) {
t_eshelby_stress_dx(
i,
j,
m,
n) =
(t_F(
k,
i) * t_D(
k,
j,
m,
l)) * t_invH(
n,
l);
for (int ii = 0; ii != 3; ++ii)
for (int jj = 0; jj != 3; ++jj)
for (int mm = 0; mm != 3; ++mm)
for (int nn = 0; nn != 3; ++nn) {
auto &
v = t_eshelby_stress_dx(ii, jj, mm, nn);
v += t_invH(nn, ii) * t_cauchy_stress(mm, jj);
}
t_eshelby_stress_dx(
i,
j,
k,
l) *= -1;
t_energy_dx(
m,
n) = t_invH(
n,
j) * t_cauchy_stress(
m,
j);
for (int mm = 0; mm != 3; ++mm)
for (int nn = 0; nn != 3; ++nn) {
auto v = t_energy_dx(mm, nn);
for (int ii = 0; ii != 3; ++ii)
t_eshelby_stress_dx(ii, ii, mm, nn) +=
v;
}
++t_D;
++t_invH;
++t_cauchy_stress;
++t_eshelby_stress_dx;
++t_F;
}
}
template <class ELEMENT>
HookeElement::OpPostProcHookeElement<ELEMENT>::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, bool is_field_disp)
dataAtPts(data_at_pts), blockSetsPtr(block_sets_ptr),
postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts), isALE(is_ale),
isFieldDisp(is_field_disp) {}
template <class ELEMENT>
int side,
EntityType type, EntitiesFieldData::EntData &data) {
if (type != MBVERTEX) {
}
auto tensor_to_tensor = [](const auto &t1, auto &t2) {
t2(0, 0) = t1(0, 0);
t2(1, 1) = t1(1, 1);
t2(2, 2) = t1(2, 2);
t2(0, 1) = t2(1, 0) = t1(1, 0);
t2(0, 2) = t2(2, 0) = t1(2, 0);
t2(1, 2) = t2(2, 1) = t1(2, 1);
};
std::array<double, 9> def_val;
def_val.fill(0);
auto make_tag = [&](auto name, auto size) {
CHKERR postProcMesh.tag_get_handle(name, size, MB_TYPE_DOUBLE, th,
MB_TAG_CREAT | MB_TAG_SPARSE,
def_val.data());
};
auto th_stress = make_tag("STRESS", 9);
auto th_psi = make_tag("ENERGY", 1);
const int nb_integration_pts = mapGaussPts.size();
auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
auto t_H = getFTensor2FromMat<3, 3>(*dataAtPts->HMat);
dataAtPts->stiffnessMat->resize(36, 1, false);
if (type == MBTRI || type == MBQUAD) {
auto &m_field = this->getPtrFE()->mField;
CHKERR m_field.get_moab().get_adjacencies(&ent, 1, 3,
false, ents,
moab::Interface::UNION);
#ifndef NDEBUG
if (ents.empty())
"Could not find a 3D element adjacent to a given face element");
#endif
ent_3d = ents.front();
}
bool found_block = false;
for (
auto &
m : (blockSetsPtr)) {
if (
m.second.tEts.find(ent_3d) !=
m.second.tEts.end()) {
const double young =
m.second.E;
const double poisson =
m.second.PoissonRatio;
const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
t_D(0, 0, 0, 0) = t_D(1, 1, 1, 1) = t_D(2, 2, 2, 2) = 1 - poisson;
t_D(0, 1, 0, 1) = t_D(0, 2, 0, 2) = t_D(1, 2, 1, 2) =
0.5 * (1 - 2 * poisson);
t_D(0, 0, 1, 1) = t_D(1, 1, 0, 0) = t_D(0, 0, 2, 2) = t_D(2, 2, 0, 0) =
t_D(1, 1, 2, 2) = t_D(2, 2, 1, 1) = poisson;
t_D(
i,
j,
k,
l) *= coefficient;
found_block = true;
break;
}
}
if (!found_block)
"Element not found in any of material blocksets");
double detH = 0.;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
if (isFieldDisp) {
t_h(0, 0) += 1;
t_h(1, 1) += 1;
t_h(2, 2) += 1;
}
if (!isALE) {
t_small_strain_symm(
i,
j) = (t_h(
i,
j) || t_h(
j,
i)) / 2.;
} else {
t_F(
i,
j) = t_h(
i,
k) * t_invH(
k,
j);
t_small_strain_symm(
i,
j) = (t_F(
i,
j) || t_F(
j,
i)) / 2.;
++t_H;
}
t_small_strain_symm(0, 0) -= 1;
t_small_strain_symm(1, 1) -= 1;
t_small_strain_symm(2, 2) -= 1;
t_stress_symm(
i,
j) = t_D(
i,
j,
k,
l) * t_small_strain_symm(
k,
l);
tensor_to_tensor(t_stress_symm, t_stress);
const double psi = 0.5 * t_stress_symm(
i,
j) * t_small_strain_symm(
i,
j);
CHKERR postProcMesh.tag_set_data(th_psi, &mapGaussPts[gg], 1, &psi);
CHKERR postProcMesh.tag_set_data(th_stress, &mapGaussPts[gg], 1,
&t_stress(0, 0));
++t_h;
}
}
template <int S>
HookeElement::OpAnalyticalInternalStrain_dx<S>::OpAnalyticalInternalStrain_dx(
const std::string row_field,
boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
StrainFunction strain_fun)
:
OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
strainFun(strain_fun) {}
template <int S>
HookeElement::OpAnalyticalInternalStrain_dx<S>::iNtegrate(
EntData &row_data) {
&
v(r + 0), &
v(r + 1), &
v(r + 2));
};
const int nb_integration_pts = getGaussPts().size2();
auto t_coords = getFTensor1CoordsAtGaussPts();
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
nF.resize(nbRows, false);
nF.clear();
const int row_nb_base_fun = row_data.
getN().size2();
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
auto t_fun_strain = strainFun(t_coords);
t_stress(
i,
j) = -t_D(
i,
j,
k,
l) * t_fun_strain(
k,
l);
auto t_nf = get_tensor1(nF, 0);
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
t_nf(
i) +=
a * t_row_diff_base(
j) * t_stress(
i,
j);
++t_row_diff_base;
++t_nf;
}
for (; rr != row_nb_base_fun; ++rr)
++t_row_diff_base;
++t_w;
++t_coords;
++t_D;
}
}
template <int S>
HookeElement::OpAnalyticalInternalAleStrain_dX<S>::
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)
:
OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
strainFun(strain_fun), matPosAtPtsPtr(mat_pos_at_pts_ptr) {}
template <int S>
MoFEMErrorCode HookeElement::OpAnalyticalInternalAleStrain_dX<S>::iNtegrate(
&
v(r + 0), &
v(r + 1), &
v(r + 2));
};
const int nb_integration_pts = getGaussPts().size2();
auto get_coords = [&]() { return getFTensor1FromMat<3>(*matPosAtPtsPtr); };
auto t_coords = get_coords();
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
nF.resize(nbRows, false);
nF.clear();
auto t_F = getFTensor2FromMat<3, 3>(*(dataAtPts->FMat));
auto &det_H = *dataAtPts->detHVec;
auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
const int row_nb_base_fun = row_data.
getN().size2();
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
auto t_fun_strain = strainFun(t_coords);
t_stress(
i,
j) = -t_D(
i,
j,
k,
l) * t_fun_strain(
k,
l);
t_eshelby_stress(
i,
j) = -t_F(
k,
i) * t_stress(
k,
j);
double a = t_w * vol * det_H[gg];
auto t_nf = get_tensor1(nF, 0);
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
t_nf(
i) +=
a * t_row_diff_base_pulled(
j) * t_eshelby_stress(
i,
j);
++t_row_diff_base;
++t_nf;
}
for (; rr != row_nb_base_fun; ++rr)
++t_row_diff_base;
++t_w;
++t_coords;
++t_F;
++t_invH;
++t_D;
}
}
template <int S>
HookeElement::OpAnalyticalInternalAleStrain_dx<S>::
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)
:
OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
strainFun(strain_fun), matPosAtPtsPtr(mat_pos_at_pts_ptr) {}
template <int S>
MoFEMErrorCode HookeElement::OpAnalyticalInternalAleStrain_dx<S>::iNtegrate(
&
v(r + 0), &
v(r + 1), &
v(r + 2));
};
const int nb_integration_pts = getGaussPts().size2();
auto get_coords = [&]() { return getFTensor1FromMat<3>(*matPosAtPtsPtr); };
auto t_coords = get_coords();
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
nF.resize(nbRows, false);
nF.clear();
auto &det_H = *dataAtPts->detHVec;
auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
const int row_nb_base_fun = row_data.
getN().size2();
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
auto t_fun_strain = strainFun(t_coords);
t_stress(
i,
j) = -t_D(
i,
j,
k,
l) * t_fun_strain(
k,
l);
double a = t_w * vol * det_H[gg];
auto t_nf = get_tensor1(nF, 0);
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
t_nf(
i) +=
a * t_row_diff_base_pulled(
j) * t_stress(
i,
j);
++t_row_diff_base;
++t_nf;
}
for (; rr != row_nb_base_fun; ++rr)
++t_row_diff_base;
++t_w;
++t_coords;
++t_invH;
++t_D;
}
}
#endif
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.
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
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...
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
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
Assemble mass matrix for elastic element TODO: CHANGE FORMULA *.