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

Operators and data structures for calculating internal stress.

Operators and data structures for calculating internal stress

/**
* \file HookeInternalStressElement.hpp
* \example HookeInternalStressElement.hpp
*
* \brief Operators and data structures for calculating internal stress
*
*/
/*
* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
#ifndef __HOOKE_INTERNAL_STRESS_ELEMENT_HPP__
#define __HOOKE_INTERNAL_STRESS_ELEMENT_HPP__
struct DataAtIntegrationPts : public HookeElement::DataAtIntegrationPts {
boost::shared_ptr<MatrixDouble> internalStressMat;
boost::shared_ptr<MatrixDouble> actualStressMat;
boost::shared_ptr<MatrixDouble> deviatoricStressMat;
boost::shared_ptr<MatrixDouble> hydrostaticStressMat;
boost::shared_ptr<MatrixDouble> spatPosMat;
boost::shared_ptr<MatrixDouble> meshNodePosMat;
internalStressMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
actualStressMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
deviatoricStressMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
boost::shared_ptr<MatrixDouble>(new MatrixDouble());
spatPosMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
meshNodePosMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
}
};
struct OpGetInternalStress : public VolUserDataOperator {
OpGetInternalStress(const std::string row_field,
const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
moab::Interface &input_mesh, char *stress_tag_name)
: VolUserDataOperator(row_field, col_field, OPROW, true),
dataAtPts(data_at_pts), inputMesh(input_mesh),
stressTagName(stress_tag_name) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
};
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
protected:
boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
moab::Interface &inputMesh;
};
struct OpInternalStrain_dx : public OpAssemble {
OpInternalStrain_dx(const std::string row_field,
boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
: OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
dataAtPts(data_at_pts){};
MoFEMErrorCode iNtegrate(EntData &row_data);
protected:
boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
};
template <int S = 0>
struct OpGetAnalyticalInternalStress : public VolUserDataOperator {
typedef boost::function<
)
>
OpGetAnalyticalInternalStress(
const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
StrainFunction strain_fun)
: VolUserDataOperator(row_field, col_field, OPROW, true),
dataAtPts(data_at_pts), strainFun(strain_fun) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
protected:
boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
};
struct OpSaveStress : public VolUserDataOperator {
boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
map<int, BlockData>
&blockSetsPtr; // FIXME: (works only with the first block)
moab::Interface &outputMesh;
bool isALE;
double scaleFactor;
bool saveMean;
OpSaveStress(const string row_field, const string col_field,
boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
map<int, BlockData> &block_sets_ptr,
moab::Interface &output_mesh, double scale_factor,
bool save_mean = false, bool is_ale = false,
bool is_field_disp = true)
: VolUserDataOperator(row_field, col_field, OPROW, true),
dataAtPts(data_at_pts), blockSetsPtr(block_sets_ptr),
outputMesh(output_mesh), scaleFactor(scale_factor), isALE(is_ale),
isFieldDisp(is_field_disp), saveMean(save_mean){};
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
};
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;
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)
: 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) {}
MoFEMErrorCode doWork(int side, EntityType type,
DataForcesAndSourcesCore::EntData &data);
};
};
int row_side, EntityType row_type, EntData &row_data) {
const int nb_integration_pts = row_data.getN().size1();
const int val_num = 9 * nb_integration_pts;
std::vector<double> def_vals(val_num, 0.0);
Tag th_internal_stress;
CHKERR inputMesh.tag_get_handle(
stressTagName, val_num, MB_TYPE_DOUBLE, th_internal_stress,
MB_TAG_CREAT | MB_TAG_SPARSE, &*def_vals.begin());
dataAtPts->internalStressMat->resize(9, nb_integration_pts, false);
CHKERR inputMesh.tag_get_data(
th_internal_stress, &ent, 1,
&*(dataAtPts->internalStressMat->data().begin()));
}
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_internal_stress =
getFTensor2FromMat<3, 3>(*dataAtPts->internalStressMat);
// get element volume
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
nF.resize(nbRows, false);
nF.clear();
// 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) {
// 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_internal_stress(i, j);
++t_row_diff_base;
++t_nf;
}
for (; rr != row_nb_base_fun; ++rr)
++t_row_diff_base;
++t_w;
++t_internal_stress;
}
}
template <int S>
int row_side, EntityType row_type, EntData &row_data) {
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);
};
const int nb_integration_pts = getGaussPts().size2();
auto t_coords = getFTensor1CoordsAtGaussPts();
// elastic stiffness tensor (4th rank tensor with minor and major
// symmetry)
MAT_TO_DDG(dataAtPts->stiffnessMat));
dataAtPts->internalStressMat->resize(9, nb_integration_pts, false);
auto t_internal_stress =
getFTensor2FromMat<3, 3>(*dataAtPts->internalStressMat);
FTensor::Tensor2_symmetric<double, 3> t_internal_stress_symm;
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
auto t_fun_strain = strainFun(t_coords);
t_internal_stress_symm(i, j) = -t_D(i, j, k, l) * t_fun_strain(k, l);
tensor_to_tensor(t_internal_stress_symm, t_internal_stress);
++t_coords;
++t_D;
++t_internal_stress;
}
}
int row_side, EntityType row_type, EntData &row_data) {
if (row_type != MBVERTEX) {
}
const EntityHandle ent = getFEEntityHandle();
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);
};
auto tensor_to_vector = [](const auto &t, auto &v) {
v(0) = t(0, 0);
v(1) = t(1, 1);
v(2) = t(2, 2);
v(3) = t(0, 1);
v(4) = t(0, 2);
v(5) = t(1, 2);
};
auto get_tag_handle = [&](auto name, auto size) {
Tag th;
std::vector<double> def_vals(size, 0.0);
CHKERR outputMesh.tag_get_handle(name, size, MB_TYPE_DOUBLE, th,
MB_TAG_CREAT | MB_TAG_SPARSE,
def_vals.data());
return th;
};
const int nb_integration_pts = row_data.getN().size1();
auto th_internal_stress =
get_tag_handle("INTERNAL_STRESS", 9 * nb_integration_pts);
auto th_actual_stress =
get_tag_handle("ACTUAL_STRESS", 9 * nb_integration_pts);
auto th_deviatoric_stress =
get_tag_handle("DEVIATORIC_STRESS", 9 * nb_integration_pts);
auto th_hydrostatic_stress =
get_tag_handle("HYDROSTATIC_STRESS", 9 * nb_integration_pts);
auto th_internal_stress_mean = get_tag_handle("MED_INTERNAL_STRESS", 9);
auto th_actual_stress_mean = get_tag_handle("MED_ACTUAL_STRESS", 9);
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));
for (auto &m : (blockSetsPtr)) {
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;
break; // FIXME: calculates only first block
}
double detH = 0.;
auto t_internal_stress =
getFTensor2FromMat<3, 3>(*dataAtPts->internalStressMat);
dataAtPts->actualStressMat->resize(9, nb_integration_pts, false);
auto t_actual_stress = getFTensor2FromMat<3, 3>(*dataAtPts->actualStressMat);
dataAtPts->deviatoricStressMat->resize(9, nb_integration_pts, false);
auto t_deviatoric_stress =
getFTensor2FromMat<3, 3>(*dataAtPts->deviatoricStressMat);
dataAtPts->hydrostaticStressMat->resize(9, nb_integration_pts, false);
auto t_hydrostatic_stress =
getFTensor2FromMat<3, 3>(*dataAtPts->hydrostaticStressMat);
FTensor::Tensor2<double, 3, 3> t_internal_stress_mean;
FTensor::Tensor2<double, 3, 3> t_actual_stress_mean;
t_internal_stress_mean(i, j) = 0.;
t_actual_stress_mean(i, j) = 0.;
auto t_w = getFTensor0IntegrationWeight();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
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;
}
if (isALE || !isFieldDisp) {
for (auto ii : {0, 1, 2}) {
t_small_strain_symm(ii, ii) -= 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);
t_actual_stress(i, j) = t_stress(i, j) + t_internal_stress(i, j);
t_actual_stress(i, j) *= scaleFactor;
t_internal_stress(i, j) *= scaleFactor;
double hydrostatic_pressure =
(t_actual_stress(0, 0) + t_actual_stress(1, 1) +
t_actual_stress(2, 2)) /
3.;
t_hydrostatic_stress(i, j) = 0.;
for (auto ii : {0, 1, 2}) {
t_hydrostatic_stress(ii, ii) += hydrostatic_pressure;
}
t_deviatoric_stress(i, j) = t_actual_stress(i, j);
for (auto ii : {0, 1, 2}) {
t_deviatoric_stress(ii, ii) -= hydrostatic_pressure;
}
t_actual_stress_mean(i, j) += t_w * t_actual_stress(i, j);
t_internal_stress_mean(i, j) += t_w * t_internal_stress(i, j);
++t_w;
++t_h;
++t_actual_stress;
++t_internal_stress;
}
if (saveMean) {
VectorDouble vec_actual_stress_mean;
vec_actual_stress_mean.resize(9, false);
vec_actual_stress_mean.clear();
VectorDouble vec_internal_stress_mean;
vec_internal_stress_mean.resize(9, false);
vec_internal_stress_mean.clear();
tensor_to_vector(t_actual_stress_mean, vec_actual_stress_mean);
tensor_to_vector(t_internal_stress_mean, vec_internal_stress_mean);
CHKERR outputMesh.tag_set_data(th_actual_stress_mean, &ent, 1,
&*(vec_actual_stress_mean.data().begin()));
CHKERR outputMesh.tag_set_data(th_internal_stress_mean, &ent, 1,
&*(vec_internal_stress_mean.data().begin()));
} else {
CHKERR outputMesh.tag_set_data(
th_internal_stress, &ent, 1,
&*(dataAtPts->internalStressMat->data().begin()));
CHKERR outputMesh.tag_set_data(
th_actual_stress, &ent, 1,
&*(dataAtPts->actualStressMat->data().begin()));
CHKERR outputMesh.tag_set_data(
th_hydrostatic_stress, &ent, 1,
&*(dataAtPts->hydrostaticStressMat->data().begin()));
CHKERR outputMesh.tag_set_data(
th_deviatoric_stress, &ent, 1,
&*(dataAtPts->deviatoricStressMat->data().begin()));
}
}
template <class ELEMENT>
int side, EntityType type, DataForcesAndSourcesCore::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);
};
auto get_tag_handle = [&](auto name, auto size) {
Tag th;
std::vector<double> def_vals(size, 0.0);
CHKERR postProcMesh.tag_get_handle(name, size, MB_TYPE_DOUBLE, th,
MB_TAG_CREAT | MB_TAG_SPARSE,
def_vals.data());
return th;
};
auto th_stress = get_tag_handle("STRESS", 9);
auto th_strain = get_tag_handle("STRAIN", 9);
auto th_psi = get_tag_handle("ENERGY", 1);
auto th_disp = get_tag_handle("DISPLACEMENT", 3);
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));
for (auto &m : (blockSetsPtr)) {
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;
break; // FIXME: calculates only first block
}
double detH = 0.;
auto t_spat_pos = getFTensor1FromMat<3>(*dataAtPts->spatPosMat);
auto t_mesh_node_pos = getFTensor1FromMat<3>(*dataAtPts->meshNodePosMat);
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);
tensor_to_tensor(t_small_strain_symm, t_small_strain);
t_disp(i) = t_spat_pos(i) - t_mesh_node_pos(i);
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));
CHKERR postProcMesh.tag_set_data(th_strain, &mapGaussPts[gg], 1,
&t_small_strain(0, 0));
CHKERR postProcMesh.tag_set_data(th_disp, &mapGaussPts[gg], 1, &t_disp(0));
++t_h;
++t_spat_pos;
++t_mesh_node_pos;
}
}
#endif
#define MAT_TO_DDG(SM)
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
#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
VolumeElementForcesAndSourcesCore::UserDataOperator VolUserDataOperator
Definition: elasticity.cpp:76
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
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
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
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
const double r
rate factor
constexpr double t
plate stiffness
Definition: plate.cpp:59
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
boost::function< FTensor::Tensor2_symmetric< double, 3 >(FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, 3 > &t_coords) > StrainFunction
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
bool isALE
bool isFieldDisp
map< int, BlockData > & blockSetsPtr
moab::Interface & postProcMesh
std::vector< EntityHandle > & mapGaussPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
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....
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
@ OPROW
operator doWork function is executed on FE rows