#ifndef __ADOLC_MATERIAL_MODELS_HPP
#define __ADOLC_MATERIAL_MODELS_HPP
#ifndef WITH_ADOL_C
#error "MoFEM need to be compiled with ADOL-C"
#endif
template <int DIM> struct J2Plasticity;
template <> struct J2Plasticity<3> : public ClosestPointProjection {
J2Plasticity() : ClosestPointProjection() {
nbInternalVariables = 7;
activeVariablesW.resize(12 + nbInternalVariables, false);
}
enum ModelParams {
LAMBDA,
MU, ISOH, KINK, PHI, SIGMA_Y, LAST_PARAM };
using ParamsArray = std::array<double, LAST_PARAM>;
boost::shared_ptr<ParamsArray> defaultParamsArrayPtr = nullptr;
boost::shared_ptr<ParamsArray> paramsArrayPtr = nullptr;
boost::shared_ptr<ParamsArray> oldParamsArrayPtr = nullptr;
PETSC_NULL);
0);
MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
defaultParamsArrayPtr = boost::make_shared<ParamsArray>();
(*defaultParamsArrayPtr)[
MU] =
mu;
(*defaultParamsArrayPtr)[
H] =
H;
(*defaultParamsArrayPtr)[
K] =
K;
(*defaultParamsArrayPtr)[PHI] =
phi;
(*defaultParamsArrayPtr)[SIGMA_Y] =
sigmaY;
MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
<<
"LAMBDA: " << (*defaultParamsArrayPtr)[
LAMBDA];
MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
<<
"MU: " << (*defaultParamsArrayPtr)[
MU];
MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
<< "H: " << (*defaultParamsArrayPtr)[ISOH];
MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
<< "K: " << (*defaultParamsArrayPtr)[KINK];
MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
<< "PHI: " << (*defaultParamsArrayPtr)[PHI];
MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
<< "SIGMA_Y: " << (*defaultParamsArrayPtr)[SIGMA_Y];
paramsArrayPtr = boost::make_shared<ParamsArray>();
std::copy(defaultParamsArrayPtr->begin(), defaultParamsArrayPtr->end(),
paramsArrayPtr->begin());
oldParamsArrayPtr = boost::make_shared<ParamsArray>();
std::fill(oldParamsArrayPtr->begin(), oldParamsArrayPtr->end(), 0);
};
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string block_name,
Sev sev) {
OpMatBlocks(boost::shared_ptr<ParamsArray> params_array_ptr,
boost::shared_ptr<ParamsArray> def_params_array_otr,
std::vector<const CubitMeshSets *> meshset_vec_ptr)
:
OP(
NOSPACE, OP::OPSPACE), paramsArrayPtr(params_array_ptr),
defParamsArray(def_params_array_otr) {
"Can not get data from block");
}
for (auto &b : blockData) {
if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
std::copy(b.paramsArray.begin(), b.paramsArray.end(),
paramsArrayPtr->begin());
}
}
std::copy(defParamsArray->begin(), defParamsArray->end(),
paramsArrayPtr->begin());
}
private:
boost::shared_ptr<ParamsArray> paramsArrayPtr;
boost::shared_ptr<ParamsArray> defParamsArray;
ParamsArray paramsArray;
};
std::vector<BlockData> blockData;
std::vector<const CubitMeshSets *> meshset_vec_ptr,
for (
auto m : meshset_vec_ptr) {
std::vector<double> block_data;
CHKERR m->getAttributes(block_data);
if (block_data.size() != 6) {
"Expected that block has two attribute");
}
auto get_block_ents = [&]() {
m_field.
get_moab().get_entities_by_handle(
m->meshset, ents,
true);
return ents;
};
blockData.push_back({*defParamsArray, get_block_ents()});
std::copy(block_data.begin(), block_data.end(),
blockData.back().paramsArray.begin());
<<
"Young modulus: " << (blockData.back().paramsArray)[
LAMBDA];
<<
"Poisson ratio: " << (blockData.back().paramsArray)[
MU];
blockData.back().paramsArray[
LAMBDA] =
<<
"LAMBDA: " << (blockData.back().paramsArray)[
LAMBDA];
<<
"MU: " << (blockData.back().paramsArray)[
MU];
<< "H: " << (blockData.back().paramsArray)[ISOH];
<< "K: " << (blockData.back().paramsArray)[KINK];
<< "PHI: " << (blockData.back().paramsArray)[PHI];
<< "SIGMA_Y: " << (blockData.back().paramsArray)[SIGMA_Y];
}
}
};
auto cubit_meshsets_vec_ptr =
m_field.
getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
std::regex((boost::format("%s(.*)") % block_name).str()));
pip.push_back(new OpMatBlocks(paramsArrayPtr, defaultParamsArrayPtr,
m_field, sev, cubit_meshsets_vec_ptr));
}
MoFEMErrorCode setParams(
short tag,
bool &recalculate_elastic_tangent) {
switch (tag) {
set_param_vec(tapesTags[tag], paramsArrayPtr->size(),
paramsArrayPtr->data());
for (auto p = 0; p != LAST_PARAM; p++) {
if ((*paramsArrayPtr)[p] != (*oldParamsArrayPtr)[p]) {
recalculate_elastic_tangent = true;
std::copy(paramsArrayPtr->begin(), paramsArrayPtr->end(),
oldParamsArrayPtr->begin());
break;
}
}
} break;
case TypesTags::Y:
break;
default:
"Unknown tag for setParams");
}
}
auto p_lambda = mkparam(
lambda);
auto p_phi = mkparam(
phi);
auto p_sigma_y = mkparam(
sigmaY);
auto t_voight_total_strain = getFTensor1FromPtr<6>(&a_sTrain[0]);
auto t_voight_plastic_strain = getFTensor1FromPtr<6>(&a_plasticStrain[0]);
auto t_voight_internal_tensor =
getFTensor1FromPtr<6>(&a_internalVariables[1]);
tTotalStrain(
i,
j) = t_voight_strain_op(
i,
j, Z) * t_voight_total_strain(Z);
t_voight_strain_op(
i,
j, Z) * t_voight_plastic_strain(Z);
t_voight_strain_op(
i,
j, Z) * t_voight_internal_tensor(Z);
tElasticStrain(
i,
j) = tTotalStrain(
i,
j) - tPlasticStrain(
i,
j);
tR = tElasticStrain(
i,
i);
a_w = 0.5 * p_lambda * tR * tR +
p_mu * (tElasticStrain(
i,
j) * tElasticStrain(
i,
j));
a_w +=
a_internalVariables[0] * p_sigma_y +
0.5 * p_phi * p_H * (a_internalVariables[0] * a_internalVariables[0]);
a_w += (0.5 * (1 - p_phi) * p_K) *
(tInternalTensor(
i,
j) * tInternalTensor(
i,
j));
}
auto t_vioght_stress = getFTensor1FromPtr<6>(&a_sTress[0]);
auto t_back_stress = getFTensor1FromPtr<6>(&a_internalFluxes[1]);
tStress(
i,
j) = t_voight_stress_op(
i,
j, Z) * t_vioght_stress(Z);
tBackStress(
i,
j) = t_voight_stress_op(
i,
j, Z) * t_back_stress(Z);
constexpr
double third = boost::math::constants::third<double>();
tDeviator(
i,
j) = tStress(
i,
j) - tR *
t_kd(
i,
j) - tBackStress(
i,
j);
constexpr
double c = 3. / 2.;
f =
c * tDeviator(
i,
j) * tDeviator(
i,
j);
}
a_y = sqrt(
f) - a_internalFluxes[0];
}
a_h = sqrt(
f) - a_internalFluxes[0];
}
};
template <> struct J2Plasticity<2> : public J2Plasticity<3> {
auto p_lambda = mkparam(
lambda);
auto p_phi = mkparam(
phi);
auto p_sigma_y = mkparam(
sigmaY);
auto t_voight_total_strain = getFTensor1FromPtr<6>(&a_sTrain[0]);
auto t_voight_plastic_strain = getFTensor1FromPtr<6>(&a_plasticStrain[0]);
auto t_voight_internal_tensor =
getFTensor1FromPtr<6>(&a_internalVariables[1]);
tTotalStrain(
i,
j) = t_voight_strain_op(
i,
j, Z) * t_voight_total_strain(Z);
t_voight_strain_op(
i,
j, Z) * t_voight_plastic_strain(Z);
t_voight_strain_op(
i,
j, Z) * t_voight_internal_tensor(Z);
tElasticStrain(
i,
j) = tTotalStrain(
i,
j) - tPlasticStrain(
i,
j);
tElasticStrain(2, 2) =
(-tElasticStrain(0, 0) * p_lambda - tElasticStrain(1, 1) * p_lambda) /
(p_lambda + 2 * p_mu);
tR = tElasticStrain(
i,
i);
a_w = 0.5 * p_lambda * tR * tR +
p_mu * (tElasticStrain(
i,
j) * tElasticStrain(
i,
j));
a_w +=
a_internalVariables[0] * p_sigma_y +
0.5 * p_phi * p_H * (a_internalVariables[0] * a_internalVariables[0]);
a_w += (0.5 * (1 - p_phi) * p_K) *
(tInternalTensor(
i,
j) * tInternalTensor(
i,
j));
}
};
struct ParaboloidalPlasticity : public ClosestPointProjection {
}
};
PETSC_NULL);
0);
MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
(*defaultParamsArrayPtr)[
MU] =
mu;
(*defaultParamsArrayPtr)[
NUP] =
nup;
(*defaultParamsArrayPtr)[
HT] =
Ht;
(*defaultParamsArrayPtr)[
HC] =
Hc;
(*defaultParamsArrayPtr)[
NT] =
nt;
(*defaultParamsArrayPtr)[
NC] =
nc;
MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
<<
"LAMBDA: " << (*defaultParamsArrayPtr)[
LAMBDA];
MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
<<
"MU: " << (*defaultParamsArrayPtr)[
MU];
MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
<<
"NUP: " << (*defaultParamsArrayPtr)[
NUP];
MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
<<
"SIGMA_YT: " << (*defaultParamsArrayPtr)[
SIGMA_YT];
MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
<<
"SIGMA_YC: " << (*defaultParamsArrayPtr)[
SIGMA_YC];
MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
<<
"HT: " << (*defaultParamsArrayPtr)[
HT];
MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
<<
"HC: " << (*defaultParamsArrayPtr)[
HC];
MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
<<
"NT: " << (*defaultParamsArrayPtr)[
NT];
MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
<<
"NC: " << (*defaultParamsArrayPtr)[
NC];
};
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string block_name,
Sev sev) {
return 0;
}
recalculate_elastic_tangent = true;
return 0;
}
auto t_voight_total_strain = getFTensor1FromPtr<6>(&
a_sTrain[0]);
t_voight_strain_op(
i,
j,
Z) * t_voight_plastic_strain(
Z);
}
auto t_vioght_stress = getFTensor1FromPtr<6>(&
a_sTress[0]);
tStress(
i,
j) = t_voight_stress_op(
i,
j,
Z) * t_vioght_stress(
Z);
}
}
double alpha =
}
};
template <typename MODEL>
std::array<int, ClosestPointProjection::LAST_TAPE> tapes_tags = {0, 1, 2}) {
auto cp_ptr = boost::make_shared<MODEL>();
cp_ptr->getDefaultMaterialParameters();
cp_ptr->createMatAVecR();
cp_ptr->snesCreate();
cp_ptr->recordTapes();
cp_ptr->tapesTags = tapes_tags;
return cp_ptr;
}
}
#endif //__ADOLC_MATERIAL_MODELS_HPP