9 #ifndef __ADOLC_MATERIAL_MODELS_HPP
10 #define __ADOLC_MATERIAL_MODELS_HPP
13 #error "MoFEM need to be compiled with ADOL-C"
76 nbInternalVariables = 7;
77 activeVariablesW.resize(12 + nbInternalVariables,
false);
104 boost::shared_ptr<ParamsArray> defaultParamsArrayPtr =
nullptr;
105 boost::shared_ptr<ParamsArray> paramsArrayPtr =
nullptr;
106 boost::shared_ptr<ParamsArray> oldParamsArrayPtr =
nullptr;
127 MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
129 MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
135 defaultParamsArrayPtr = boost::make_shared<ParamsArray>();
137 (*defaultParamsArrayPtr)[
MU] =
mu;
138 (*defaultParamsArrayPtr)[
H] =
H;
139 (*defaultParamsArrayPtr)[
K] =
K;
140 (*defaultParamsArrayPtr)[PHI] =
phi;
141 (*defaultParamsArrayPtr)[SIGMA_Y] =
sigmaY;
143 MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
144 <<
"LAMBDA: " << (*defaultParamsArrayPtr)[
LAMBDA];
145 MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
146 <<
"MU: " << (*defaultParamsArrayPtr)[
MU];
147 MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
148 <<
"H: " << (*defaultParamsArrayPtr)[ISOH];
149 MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
150 <<
"K: " << (*defaultParamsArrayPtr)[KINK];
151 MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
152 <<
"PHI: " << (*defaultParamsArrayPtr)[PHI];
153 MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
154 <<
"SIGMA_Y: " << (*defaultParamsArrayPtr)[SIGMA_Y];
156 paramsArrayPtr = boost::make_shared<ParamsArray>();
157 std::copy(defaultParamsArrayPtr->begin(), defaultParamsArrayPtr->end(),
158 paramsArrayPtr->begin());
159 oldParamsArrayPtr = boost::make_shared<ParamsArray>();
160 std::fill(oldParamsArrayPtr->begin(), oldParamsArrayPtr->end(), 0);
176 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
177 std::string block_name,
Sev sev) {
184 OpMatBlocks(boost::shared_ptr<ParamsArray> params_array_ptr,
185 boost::shared_ptr<ParamsArray> def_params_array_otr,
187 std::vector<const CubitMeshSets *> meshset_vec_ptr)
188 :
OP(
NOSPACE, OP::OPSPACE), paramsArrayPtr(params_array_ptr),
189 defParamsArray(def_params_array_otr) {
191 "Can not get data from block");
198 for (
auto &b : blockData) {
200 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
201 std::copy(b.paramsArray.begin(), b.paramsArray.end(),
202 paramsArrayPtr->begin());
207 std::copy(defParamsArray->begin(), defParamsArray->end(),
208 paramsArrayPtr->begin());
214 boost::shared_ptr<ParamsArray> paramsArrayPtr;
215 boost::shared_ptr<ParamsArray> defParamsArray;
222 std::vector<BlockData> blockData;
226 std::vector<const CubitMeshSets *> meshset_vec_ptr,
230 for (
auto m : meshset_vec_ptr) {
232 std::vector<double> block_data;
233 CHKERR m->getAttributes(block_data);
234 if (block_data.size() != 6) {
236 "Expected that block has two attribute");
238 auto get_block_ents = [&]() {
241 m_field.
get_moab().get_entities_by_handle(
m->meshset, ents,
true);
245 blockData.push_back({*defParamsArray, get_block_ents()});
246 std::copy(block_data.begin(), block_data.end(),
247 blockData.back().paramsArray.begin());
255 <<
"Young modulus: " << (blockData.back().paramsArray)[
LAMBDA];
257 <<
"Poisson ratio: " << (blockData.back().paramsArray)[
MU];
259 blockData.back().paramsArray[
LAMBDA] =
264 <<
"LAMBDA: " << (blockData.back().paramsArray)[
LAMBDA];
266 <<
"MU: " << (blockData.back().paramsArray)[
MU];
268 <<
"H: " << (blockData.back().paramsArray)[ISOH];
270 <<
"K: " << (blockData.back().paramsArray)[KINK];
272 <<
"PHI: " << (blockData.back().paramsArray)[PHI];
274 <<
"SIGMA_Y: " << (blockData.back().paramsArray)[SIGMA_Y];
280 auto cubit_meshsets_vec_ptr =
281 m_field.
getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
282 std::regex((boost::format(
"%s(.*)") % block_name).str()));
284 pip.push_back(
new OpMatBlocks(paramsArrayPtr, defaultParamsArrayPtr,
285 m_field, sev, cubit_meshsets_vec_ptr));
301 set_param_vec(tapesTags[tag], paramsArrayPtr->size(),
302 paramsArrayPtr->data());
303 for (
auto p = 0; p != LAST_PARAM; p++) {
304 if ((*paramsArrayPtr)[p] != (*oldParamsArrayPtr)[p]) {
305 recalculate_elastic_tangent =
true;
306 std::copy(paramsArrayPtr->begin(), paramsArrayPtr->end(),
307 oldParamsArrayPtr->begin());
318 "Unknown tag for setParams");
328 auto p_lambda = mkparam(
lambda);
329 auto p_mu = mkparam(
mu);
330 auto p_H = mkparam(
H);
331 auto p_K = mkparam(
K);
332 auto p_phi = mkparam(
phi);
333 auto p_sigma_y = mkparam(
sigmaY);
339 auto t_voight_total_strain = getFTensor1FromPtr<6>(&a_sTrain[0]);
340 auto t_voight_plastic_strain = getFTensor1FromPtr<6>(&a_plasticStrain[0]);
341 auto t_voight_internal_tensor =
342 getFTensor1FromPtr<6>(&a_internalVariables[1]);
345 tTotalStrain(
i,
j) = t_voight_strain_op(
i,
j, Z) * t_voight_total_strain(Z);
346 tPlasticStrain(
i,
j) =
347 t_voight_strain_op(
i,
j, Z) * t_voight_plastic_strain(Z);
348 tInternalTensor(
i,
j) =
349 t_voight_strain_op(
i,
j, Z) * t_voight_internal_tensor(Z);
352 tElasticStrain(
i,
j) = tTotalStrain(
i,
j) - tPlasticStrain(
i,
j);
353 tR = tElasticStrain(
i,
i);
356 a_w = 0.5 * p_lambda * tR * tR +
357 p_mu * (tElasticStrain(
i,
j) * tElasticStrain(
i,
j));
360 a_internalVariables[0] * p_sigma_y +
361 0.5 * p_phi * p_H * (a_internalVariables[0] * a_internalVariables[0]);
363 a_w += (0.5 * (1 - p_phi) * p_K) *
364 (tInternalTensor(
i,
j) * tInternalTensor(
i,
j));
376 auto t_vioght_stress = getFTensor1FromPtr<6>(&a_sTress[0]);
377 auto t_back_stress = getFTensor1FromPtr<6>(&a_internalFluxes[1]);
378 tStress(
i,
j) = t_voight_stress_op(
i,
j, Z) * t_vioght_stress(Z);
379 tBackStress(
i,
j) = t_voight_stress_op(
i,
j, Z) * t_back_stress(Z);
381 constexpr
double third = boost::math::constants::third<double>();
383 tR = tStress(
i,
i) / 3.;
384 tDeviator(
i,
j) = tStress(
i,
j) - tR *
t_kd(
i,
j) - tBackStress(
i,
j);
387 constexpr
double c = 3. / 2.;
388 f =
c * tDeviator(
i,
j) * tDeviator(
i,
j);
396 a_y = sqrt(
f) - a_internalFluxes[0];
404 a_h = sqrt(
f) - a_internalFluxes[0];
416 auto p_lambda = mkparam(
lambda);
417 auto p_mu = mkparam(
mu);
418 auto p_H = mkparam(
H);
419 auto p_K = mkparam(
K);
420 auto p_phi = mkparam(
phi);
421 auto p_sigma_y = mkparam(
sigmaY);
424 auto t_voight_total_strain = getFTensor1FromPtr<6>(&a_sTrain[0]);
425 auto t_voight_plastic_strain = getFTensor1FromPtr<6>(&a_plasticStrain[0]);
426 auto t_voight_internal_tensor =
427 getFTensor1FromPtr<6>(&a_internalVariables[1]);
429 tTotalStrain(
i,
j) = t_voight_strain_op(
i,
j, Z) * t_voight_total_strain(Z);
430 tPlasticStrain(
i,
j) =
431 t_voight_strain_op(
i,
j, Z) * t_voight_plastic_strain(Z);
432 tInternalTensor(
i,
j) =
433 t_voight_strain_op(
i,
j, Z) * t_voight_internal_tensor(Z);
435 tElasticStrain(
i,
j) = tTotalStrain(
i,
j) - tPlasticStrain(
i,
j);
439 tElasticStrain(2, 2) =
440 (-tElasticStrain(0, 0) * p_lambda - tElasticStrain(1, 1) * p_lambda) /
441 (p_lambda + 2 * p_mu);
444 tR = tElasticStrain(
i,
i);
446 a_w = 0.5 * p_lambda * tR * tR +
447 p_mu * (tElasticStrain(
i,
j) * tElasticStrain(
i,
j));
452 a_internalVariables[0] * p_sigma_y +
453 0.5 * p_phi * p_H * (a_internalVariables[0] * a_internalVariables[0]);
455 a_w += (0.5 * (1 - p_phi) * p_K) *
456 (tInternalTensor(
i,
j) * tInternalTensor(
i,
j));
599 MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
601 MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
609 (*defaultParamsArrayPtr)[
MU] =
mu;
610 (*defaultParamsArrayPtr)[
NUP] =
nup;
613 (*defaultParamsArrayPtr)[
HT] =
Ht;
614 (*defaultParamsArrayPtr)[
HC] =
Hc;
615 (*defaultParamsArrayPtr)[
NT] =
nt;
616 (*defaultParamsArrayPtr)[
NC] =
nc;
618 MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
619 <<
"LAMBDA: " << (*defaultParamsArrayPtr)[
LAMBDA];
620 MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
621 <<
"MU: " << (*defaultParamsArrayPtr)[
MU];
622 MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
623 <<
"NUP: " << (*defaultParamsArrayPtr)[
NUP];
624 MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
625 <<
"SIGMA_YT: " << (*defaultParamsArrayPtr)[
SIGMA_YT];
626 MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
627 <<
"SIGMA_YC: " << (*defaultParamsArrayPtr)[
SIGMA_YC];
628 MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
629 <<
"HT: " << (*defaultParamsArrayPtr)[
HT];
630 MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
631 <<
"HC: " << (*defaultParamsArrayPtr)[
HC];
632 MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
633 <<
"NT: " << (*defaultParamsArrayPtr)[
NT];
634 MOFEM_LOG(
"ADOLCPlasticityWold", Sev::inform)
635 <<
"NC: " << (*defaultParamsArrayPtr)[
NC];
642 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
643 std::string block_name,
Sev sev) {
648 recalculate_elastic_tangent =
true;
657 auto t_voight_total_strain = getFTensor1FromPtr<6>(&
a_sTrain[0]);
658 auto t_voight_plastic_strain = getFTensor1FromPtr<6>(&
a_plasticStrain[0]);
662 t_voight_strain_op(
i,
j,
Z) * t_voight_plastic_strain(
Z);
684 auto t_vioght_stress = getFTensor1FromPtr<6>(&
a_sTress[0]);
685 tStress(
i,
j) = t_voight_stress_op(
i,
j,
Z) * t_vioght_stress(
Z);
720 template <
typename MODEL>
722 std::array<int, ClosestPointProjection::LAST_TAPE> tapes_tags = {0, 1, 2}) {
723 auto cp_ptr = boost::make_shared<MODEL>();
724 cp_ptr->getDefaultMaterialParameters();
725 cp_ptr->createMatAVecR();
726 cp_ptr->snesCreate();
727 cp_ptr->recordTapes();
728 cp_ptr->tapesTags = tapes_tags;
734 #endif //__ADOLC_MATERIAL_MODELS_HPP