v0.13.1
Loading...
Searching...
No Matches
HookeElement.hpp

Operators and data structures for linear elastic analysis.

Operators and data structures for linear elastic analysis

Implemention of operators for Hooke material. Implementation is extended to the case when the mesh is moving as results of topological changes, also the calculation of material forces and associated tangent matrices are added to implementation.

In other words, spatial deformation is small but topological changes large.

/**
* \file HookeElement.hpp
* \example HookeElement.hpp
*
* \brief Operators and data structures for linear elastic
* analysis
*
* Implemention of operators for Hooke material. Implementation is extended to
* the case when the mesh is moving as results of topological changes, also the
* calculation of material forces and associated tangent matrices are added to
* implementation.
*
* In other words, spatial deformation is small but topological changes large.
*/
#ifndef __HOOKE_ELEMENT_HPP
#define __HOOKE_ELEMENT_HPP
#ifndef __BASICFINITEELEMENTS_HPP__
#endif // __BASICFINITEELEMENTS_HPP__
#ifndef __NONLINEAR_ELASTIC_HPP
/** \brief data for calculation heat conductivity and heat capacity elements
* \ingroup nonlinear_elastic_elem
*/
struct BlockData {
int iD;
double E;
double PoissonRatio;
Range tEts; ///< constrains elements in block set
};
};
#endif // __NONLINEAR_ELASTIC_HPP
/** \brief structure grouping operators and data used for calculation of
* nonlinear elastic element \ingroup nonlinear_elastic_elem
*
* In order to assemble matrices and right hand vectors, the loops over
* elements, entities over that elements and finally loop over integration
* points are executed.
*
* Following implementation separate those three categories of loops and to each
* loop attach operator.
*
*/
#ifndef __CONVECTIVE_MASS_ELEMENT_HPP
/** \brief data for calculation inertia forces
* \ingroup user_modules
*/
struct BlockData {
double rho0; ///< reference density
VectorDouble a0; ///< constant acceleration
Range tEts; ///< elements in block set
};
}
#endif //__CONVECTIVE_MASS_ELEMENT_HPP
struct HookeElement {
using EntData = EntitiesFieldData::EntData;
using UserDataOperator = ForcesAndSourcesCore::UserDataOperator;
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());
hMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
FMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
HMat = 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>
OpCalculateStrain(const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
private:
boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
};
OpCalculateStrainAle(const std::string row_field,
const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
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)
template <int S = 0> struct OpCalculateStress : public VolUserDataOperator {
OpCalculateStress(const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
protected:
boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
};
OpCalculateEnergy(const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
SmartPetscObj<Vec> ghost_vec = SmartPetscObj<Vec>());
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
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);
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
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);
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
protected:
boost::shared_ptr<map<int, BlockData>>
blockSetsPtr; ///< Structure keeping data about problem, like
///< material parameters
boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
};
/** * @brief Assemble mass matrix for elastic element TODO: CHANGE FORMULA *
* \f[
* {\bf{M}} = \int\limits_\Omega
* \f]
*
*/
: public VolumeElementForcesAndSourcesCore::UserDataOperator {
MatrixDouble translocK;
BlockData &dAta;
MassBlockData &massData;
boost::shared_ptr<DataAtIntegrationPts> commonData;
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)
: VolumeElementForcesAndSourcesCore::UserDataOperator(
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,
EntityType col_type,
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
&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);
// get number of integration points
// Set size can clear local tangent matrix
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;
// get integration weights
auto t_w = getFTensor0IntegrationWeight();
// integrate local matrix for entity block
for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
auto t_row_base_func = row_data.getFTensor0N(gg, 0);
// Get volume and integration weight
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;
// Next base function for column
++t_col_base_func;
}
// Next base function for row
++t_row_base_func;
}
// Next integration point for getting weight
++t_w;
}
CHKERR MatSetValues(getKSPB(), row_data, col_data, &locK(0, 0),
ADD_VALUES);
// is symmetric
if (row_type != col_type || row_side != col_side) {
translocK.resize(col_nb_dofs, row_nb_dofs, false);
noalias(translocK) = trans(locK);
CHKERR MatSetValues(getKSPB(), col_data, row_data, &translocK(0, 0),
ADD_VALUES);
}
}
};
protected:
boost::shared_ptr<map<int, BlockData>>
blockSetsPtr; ///< Structure keeping data about problem, like
///< material parameters
boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
const double rhoN; ///< exponent n in E(p) = E * (p / p_0)^n
const double rHo0; ///< p_0 reference density in E(p) = E * (p / p_0)^n
// // where p is density, E - youngs modulus
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);
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
};
struct OpAssemble : public VolUserDataOperator {
OpAssemble(const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
const char type, bool symm = false);
/**
* \brief Do calculations for give operator
* @param row_side row side number (local number) of entity on element
* @param col_side column side number (local number) of entity on element
* @param row_type type of row entity MBVERTEX, MBEDGE, MBTRI or MBTET
* @param col_type type of column entity MBVERTEX, MBEDGE, MBTRI or MBTET
* @param row_data data for row
* @param col_data data for column
* @return error code
*/
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type, EntData &row_data,
EntData &col_data);
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
protected:
// Finite element stiffness sub-matrix K_ij
MatrixDouble transK;
boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
VectorInt rowIndices;
VectorInt colIndices;
int nbRows; ///< number of dofs on rows
int nbCols; ///< number if dof on column
int nbIntegrationPts; ///< number of integration points
bool isDiag; ///< true if this block is on diagonal
virtual MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
virtual MoFEMErrorCode iNtegrate(EntData &row_data);
/**
* \brief Assemble local entity block matrix
* @param row_data row data (consist base functions on row entity)
* @param col_data column data (consist base functions on column
* entity)
* @return error code
*/
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data);
/**
* \brief Assemble local entity right-hand vector
* @param row_data row data (consist base functions on row entity)
* @param col_data column data (consist base functions on column
* entity)
* @return error code
*/
MoFEMErrorCode aSsemble(EntData &row_data);
};
struct OpRhs_dx : public OpAssemble {
OpRhs_dx(const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
protected:
MoFEMErrorCode iNtegrate(EntData &row_data);
};
template <int S = 0> struct OpLhs_dx_dx : public OpAssemble {
OpLhs_dx_dx(const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
protected:
/**
* \brief Integrate B^T D B operator
* @param row_data row data (consist base functions on row entity)
* @param col_data column data (consist base functions on column entity)
* @return error code
*/
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
};
struct OpAleRhs_dx : public OpAssemble {
OpAleRhs_dx(const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
protected:
MoFEMErrorCode iNtegrate(EntData &row_data);
};
template <int S = 0> struct OpAleLhs_dx_dx : public OpAssemble {
OpAleLhs_dx_dx(const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
protected:
/**
* \brief Integrate B^T D B operator
* @param row_data row data (consist base functions on row entity)
* @param col_data column data (consist base functions on column entity)
* @return error code
*/
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
};
template <int S = 0> struct OpAleLhs_dx_dX : public OpAssemble {
OpAleLhs_dx_dX(const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
protected:
/**
* \brief Integrate tangent stiffness for spatial momentum
* @param row_data row data (consist base functions on row entity)
* @param col_data column data (consist base functions on column entity)
* @return error code
*/
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
};
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:
/**
* \brief Integrate tangent stiffness for spatial momentum
* @param row_data row data (consist base functions on row entity)
* @param col_data column data (consist base functions on column entity)
* @return error code
*/
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
};
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:
/**
* \brief Integrate tangent stiffness for material momentum
* @param row_data row data (consist base functions on row entity)
* @param col_data column data (consist base functions on column entity)
* @return error code
*/
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
};
struct OpAleRhs_dX : public OpAssemble {
OpAleRhs_dX(const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
protected:
MoFEMErrorCode iNtegrate(EntData &row_data);
};
template <int S = 0> struct OpAleLhs_dX_dX : public OpAssemble {
OpAleLhs_dX_dX(const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
protected:
/**
* \brief Integrate tangent stiffness for material momentum
* @param row_data row data (consist base functions on row entity)
* @param col_data column data (consist base functions on column entity)
* @return error code
*/
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
};
template <int S = 0> struct OpAleLhsPre_dX_dx : public VolUserDataOperator {
OpAleLhsPre_dX_dx(const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
private:
boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
};
struct OpAleLhs_dX_dx : public OpAssemble {
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:
/**
* \brief Integrate tangent stiffness for material momentum
* @param row_data row data (consist base functions on row entity)
* @param col_data column data (consist base functions on column entity)
* @return error code
*/
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
};
template <int S> struct OpAnalyticalInternalStrain_dx : public OpAssemble {
typedef boost::function<
)
>
StrainFunction;
const std::string row_field,
boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
StrainFunction strain_fun);
protected:
MoFEMErrorCode iNtegrate(EntData &row_data);
StrainFunction strainFun;
};
template <int S> struct OpAnalyticalInternalAleStrain_dX : public OpAssemble {
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:
MoFEMErrorCode iNtegrate(EntData &row_data);
StrainFunction strainFun;
boost::shared_ptr<MatrixDouble> matPosAtPtsPtr;
};
template <int S> struct OpAnalyticalInternalAleStrain_dx : public OpAssemble {
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:
MoFEMErrorCode iNtegrate(EntData &row_data);
StrainFunction strainFun;
boost::shared_ptr<MatrixDouble> matPosAtPtsPtr;
};
template <class ELEMENT>
struct OpPostProcHookeElement : public ELEMENT::UserDataOperator {
boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
map<int, BlockData>
&blockSetsPtr; // FIXME: (works only with the first block)
moab::Interface &postProcMesh;
std::vector<EntityHandle> &mapGaussPts;
bool isALE;
bool isFieldDisp;
OpPostProcHookeElement(const string row_field,
boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
map<int, BlockData> &block_sets_ptr,
moab::Interface &post_proc_mesh,
std::vector<EntityHandle> &map_gauss_pts,
bool is_ale = false, bool is_field_disp = true);
MoFEMErrorCode doWork(int side, EntityType type,
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,
const EntityType type = MBTET,
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)
: VolUserDataOperator(row_field, col_field, OPROW, true),
dataAtPts(data_at_pts) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
template <bool D>
MoFEMErrorCode HookeElement::OpCalculateStrain<D>::doWork(int row_side,
EntityType row_type,
EntData &row_data) {
// get number of integration points
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.;
// If displacement field, not field o spatial positons is given
if (!D) {
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>
MoFEMErrorCode HookeElement::OpAleLhs_dx_dx<S>::iNtegrate(EntData &row_data,
EntData &col_data) {
// get sub-block (3x3) of local stiffens matrix, here represented by
// second order tensor
auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
&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),
&m(r + 2, c + 2));
};
// get element volume
double vol = getVolume();
// get intergrayion weights
auto t_w = getFTensor0IntegrationWeight();
// get derivatives of base functions on rows
auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
const int row_nb_base_fun = row_data.getN().size2();
// Elastic stiffness tensor (4th rank tensor with minor and major
// symmetry)
MAT_TO_DDG(dataAtPts->stiffnessMat));
auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
auto &det_H = *dataAtPts->detHVec;
// iterate over integration points
for (int gg = 0; gg != nbIntegrationPts; ++gg) {
// calculate scalar weight times element volume
double a = t_w * vol * det_H[gg];
// iterate over row base functions
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
// get sub matrix for the row
auto t_m = get_tensor2(K, 3 * rr, 0);
FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
// I mix up the indices here so that it behaves like a
// Dg. That way I don't have to have a separate wrapper
// class Christof_Expr, which simplifies things.
t_rowD(l, j, k) = t_D(i, j, k, l) * (a * t_row_diff_base_pulled(i));
// get derivatives of base functions for columns
auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
// iterate column base functions
for (int cc = 0; cc != nbCols / 3; ++cc) {
FTensor::Tensor1<double, 3> t_col_diff_base_pulled;
t_col_diff_base_pulled(j) = t_col_diff_base(i) * t_invH(i, j);
// integrate block local stiffens matrix
t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base_pulled(k);
// move to next column base function
++t_col_diff_base;
// move to next block of local stiffens matrix
++t_m;
}
// move to next row base function
++t_row_diff_base;
}
for (; rr != row_nb_base_fun; ++rr)
++t_row_diff_base;
// move to next integration weight
++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)
: VolUserDataOperator(row_field, col_field, OPROW, true),
blockSetsPtr(block_sets_ptr), dataAtPts(data_at_pts) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
template <int S>
MoFEMErrorCode HookeElement::OpCalculateHomogeneousStiffness<S>::doWork(
int row_side, EntityType row_type, EntData &row_data) {
for (auto &m : (*blockSetsPtr)) {
if (m.second.tEts.find(getFEEntityHandle()) != m.second.tEts.end()) {
dataAtPts->stiffnessMat->resize(36, 1, false);
MAT_TO_DDG(dataAtPts->stiffnessMat));
const double young = m.second.E;
const double poisson = m.second.PoissonRatio;
// coefficient used in intermediate calculation
const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
t_D(i, j, k, l) = 0.;
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)
: VolUserDataOperator(row_field, col_field, OPROW, true),
dataAtPts(data_at_pts) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
template <int S>
MoFEMErrorCode HookeElement::OpCalculateStress<S>::doWork(int row_side,
EntityType row_type,
EntData &row_data) {
// get number of integration points
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));
// elastic stiffness tensor (4th rank tensor with minor and major
// symmetry)
MAT_TO_DDG(dataAtPts->stiffnessMat));
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>
MoFEMErrorCode HookeElement::OpLhs_dx_dx<S>::iNtegrate(EntData &row_data,
EntData &col_data) {
// get sub-block (3x3) of local stiffens matrix, here represented by
// second order tensor
auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
&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),
&m(r + 2, c + 2));
};
// get element volume
double vol = getVolume();
// get intergrayion weights
auto t_w = getFTensor0IntegrationWeight();
// get derivatives of base functions on rows
auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
const int row_nb_base_fun = row_data.getN().size2();
// Elastic stiffness tensor (4th rank tensor with minor and major
// symmetry)
MAT_TO_DDG(dataAtPts->stiffnessMat));
// iterate over integration points
for (int gg = 0; gg != nbIntegrationPts; ++gg) {
// calculate scalar weight times element volume
double a = t_w * vol;
// iterate over row base functions
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
// get sub matrix for the row
auto t_m = get_tensor2(K, 3 * rr, 0);
// get derivatives of base functions for columns
auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
// I mix up the indices here so that it behaves like a
// Dg. That way I don't have to have a separate wrapper
// class Christof_Expr, which simplifies things.
t_rowD(l, j, k) = t_D(i, j, k, l) * (a * t_row_diff_base(i));
// iterate column base functions
for (int cc = 0; cc != nbCols / 3; ++cc) {
// integrate block local stiffens matrix
t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base(k);
// move to next column base function
++t_col_diff_base;
// move to next block of local stiffens matrix
++t_m;
}
// move to next row base function
++t_row_diff_base;
}
for (; rr != row_nb_base_fun; ++rr)
++t_row_diff_base;
// move to next integration weight
++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>
MoFEMErrorCode HookeElement::OpAleLhs_dx_dX<S>::iNtegrate(EntData &row_data,
EntData &col_data) {
// get sub-block (3x3) of local stiffens matrix, here represented by
// second order tensor
auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
&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),
&m(r + 2, c + 2));
};
// get element volume
double vol = getVolume();
// get intergrayion weights
auto t_w = getFTensor0IntegrationWeight();
// get derivatives of base functions on rows
auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
const int row_nb_base_fun = row_data.getN().size2();
// Elastic stiffness tensor (4th rank tensor with minor and major
// symmetry)
MAT_TO_DDG(dataAtPts->stiffnessMat));
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;
// iterate over integration points
for (int gg = 0; gg != nbIntegrationPts; ++gg) {
// calculate scalar weight times element volume
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);
// iterate over row base functions
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
// get sub matrix for the row
auto t_m = get_tensor2(K, 3 * rr, 0);
FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
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);
FTensor::Tensor3<double, 3, 3, 3> t_row_diff_base_pulled_dX;
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);
// FIXME: This operator is not implemented, doing operation by hand
// t_row_stress_dX(i, m, n) = t_row_D(i, k, l) * t_F_dX(k, l, m, n);
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);
}
// get derivatives of base functions for columns
auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
// iterate column base functions
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);
// move to next column base function
++t_col_diff_base;
// move to next block of local stiffens matrix
++t_m;
}
// move to next row base function
++t_row_diff_base;
}
for (; rr != row_nb_base_fun; ++rr)
++t_row_diff_base;
// move to next integration weight
++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>
MoFEMErrorCode HookeElement::OpAleLhs_dX_dX<S>::iNtegrate(EntData &row_data,
EntData &col_data) {
// get sub-block (3x3) of local stiffens matrix, here represented by
// second order tensor
auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
&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),
&m(r + 2, c + 2));
};
// get element volume
double vol = getVolume();
// get intergrayion weights
auto t_w = getFTensor0IntegrationWeight();
// get derivatives of base functions on rows
auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
const int row_nb_base_fun = row_data.getN().size2();
// Elastic stiffness tensor (4th rank tensor with minor and major
// symmetry)
MAT_TO_DDG(dataAtPts->stiffnessMat));
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;
// iterate over integration points
for (int gg = 0; gg != nbIntegrationPts; ++gg) {
// calculate scalar weight times element volume
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_energy_dX(k, l) +=
(t_strain(m, n) * t_D(m, n, i, j)) * t_F_dX(i, j, k, l);
t_energy_dX(k, l) /= 2.;
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;
}
// iterate over row base functions
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
// get sub matrix for the row
auto t_m = get_tensor2(K, 3 * rr, 0);
FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
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);
FTensor::Tensor3<double, 3, 3, 3> t_row_diff_base_pulled_dX;
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);
// get derivatives of base functions for columns
auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
// iterate column base functions
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);
// move to next column base function
++t_col_diff_base;
// move to next block of local stiffens matrix
++t_m;
}
// move to next row base function
++t_row_diff_base;
}
for (; rr != row_nb_base_fun; ++rr)
++t_row_diff_base;
// move to next integration weight
++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)
: VolUserDataOperator(row_field, col_field, OPROW, true),
dataAtPts(data_at_pts) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
template <int S>
MoFEMErrorCode HookeElement::OpAleLhsPre_dX_dx<S>::doWork(int row_side,
EntityType row_type,
EntData &row_data) {
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();
// Elastic stiffness tensor (4th rank tensor with minor and major
// symmetry)
MAT_TO_DDG(dataAtPts->stiffnessMat));
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)
: ELEMENT::UserDataOperator(row_field, UserDataOperator::OPROW),
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>
MoFEMErrorCode HookeElement::OpPostProcHookeElement<ELEMENT>::doWork(
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) {
Tag th;
CHKERR postProcMesh.tag_get_handle(name, size, MB_TYPE_DOUBLE, th,
MB_TAG_CREAT | MB_TAG_SPARSE,
def_val.data());
return th;
};
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);
MAT_TO_DDG(dataAtPts->stiffnessMat));
EntityHandle ent = this->getFEEntityHandle();
auto type = type_from_handle(ent);
EntityHandle ent_3d = ent;
if (type == MBTRI || type == MBQUAD) {
Range ents;
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())
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"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(i, j, k, l) = 0.;
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)
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"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 {
CHKERR invertTensor3by3(t_H, detH, t_invH);
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;
// symmetric tensors need improvement
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) {
auto get_tensor1 = [](VectorDouble &v, const int r) {
&v(r + 0), &v(r + 1), &v(r + 2));
};
const int nb_integration_pts = getGaussPts().size2();
auto t_coords = getFTensor1CoordsAtGaussPts();
// get element volume
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
nF.resize(nbRows, false);
nF.clear();
// elastic stiffness tensor (4th rank tensor with minor and major
// symmetry)
MAT_TO_DDG(dataAtPts->stiffnessMat));
// get derivatives of base functions on rows
auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
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);
// calculate scalar weight times element volume
double a = t_w * vol;
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(
EntData &row_data) {
auto get_tensor1 = [](VectorDouble &v, const int r) {
&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();
// get element volume
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
nF.resize(nbRows, false);
nF.clear();
// elastic stiffness tensor (4th rank tensor with minor and major
// symmetry)
MAT_TO_DDG(dataAtPts->stiffnessMat));
auto t_F = getFTensor2FromMat<3, 3>(*(dataAtPts->FMat));
auto &det_H = *dataAtPts->detHVec;
auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
// get derivatives of base functions on rows
auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
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);
// calculate scalar weight times element volume
double a = t_w * vol * det_H[gg];
auto t_nf = get_tensor1(nF, 0);
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
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(
EntData &row_data) {
auto get_tensor1 = [](VectorDouble &v, const int r) {
&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();
// get element volume
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
nF.resize(nbRows, false);
nF.clear();
// elastic stiffness tensor (4th rank tensor with minor and major
// symmetry)
MAT_TO_DDG(dataAtPts->stiffnessMat));
auto &det_H = *dataAtPts->detHVec;
auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
// get derivatives of base functions on rows
auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
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);
// calculate scalar weight times element volume
double a = t_w * vol * det_H[gg];
auto t_nf = get_tensor1(nF, 0);
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
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 // __HOOKE_ELEMENT_HPP
static Index< 'K', 3 > K
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)
#define MAT_TO_DDG(SM)
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData > > &block_sets_ptr)
MatrixDouble invJac
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
constexpr double a
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
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)
double D
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.
Definition: Exceptions.hpp:56
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
UBlasVector< int > VectorInt
Definition: Types.hpp:67
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1592
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.
Definition: Templates.hpp:1323
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
Definition: Templates.hpp:1342
double w(const double g, const double t)
const double r
rate factor
data for calculation inertia forces
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 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 *.