#ifndef __HOOKE_INTERNAL_STRESS_ELEMENT_HPP__
#define __HOOKE_INTERNAL_STRESS_ELEMENT_HPP__
}
};
const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
};
protected:
boost::shared_ptr<DataAtIntegrationPts>
dataAtPts;
};
boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
protected:
boost::shared_ptr<DataAtIntegrationPts>
dataAtPts;
};
template <int S = 0>
typedef boost::function<
)
>
const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
}
protected:
boost::shared_ptr<DataAtIntegrationPts>
dataAtPts;
};
boost::shared_ptr<DataAtIntegrationPts>
dataAtPts;
map<int, BlockData>
OpSaveStress(
const string row_field,
const string col_field,
boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
map<int, BlockData> &block_sets_ptr,
bool save_mean = false, bool is_ale = false,
bool is_field_disp = true)
};
template <class ELEMENT>
boost::shared_ptr<DataAtIntegrationPts>
dataAtPts;
map<int, BlockData>
boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
map<int, BlockData> &block_sets_ptr,
std::vector<EntityHandle> &map_gauss_pts,
bool is_ale = false, bool is_field_disp = true)
};
};
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;
MB_TAG_CREAT | MB_TAG_SPARSE, &*def_vals.begin());
dataAtPts->internalStressMat->resize(9, nb_integration_pts,
false);
th_internal_stress, &ent, 1,
&*(
dataAtPts->internalStressMat->data().begin()));
}
&
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);
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_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();
dataAtPts->internalStressMat->resize(9, nb_integration_pts, false);
auto t_internal_stress =
getFTensor2FromMat<3, 3>(*dataAtPts->internalStressMat);
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) {
}
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) {
};
auto get_tag_handle = [&](auto name, auto size) {
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());
};
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);
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(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;
}
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);
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 {
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.;
}
}
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) {
vec_actual_stress_mean.resize(9, false);
vec_actual_stress_mean.clear();
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>
}
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) {
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());
};
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);
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(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;
}
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 {
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);
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