1009 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1010 std::string block_name, Sev sev)
override {
1013 struct OpMatBlocks :
public ForcesAndSourcesCore::UserDataOperator {
1015 using OP = ForcesAndSourcesCore::UserDataOperator;
1017 double poissonRatio = 0.25;
1019 Tag thCompressiveYieldStress;
1020 Tag thTensionYieldStress;
1021 double youngsModulus;
1022 double yieldStrengthC;
1023 double yieldStrengthT;
1025 OpMatBlocks(boost::shared_ptr<ParamsArray> params_array_ptr,
1026 boost::shared_ptr<ParamsArray> def_params_array_ptr,
1028 std::vector<const CubitMeshSets *> meshset_vec_ptr)
1029 :
OP(
NOSPACE, OP::OPSPACE), paramsArrayPtr(params_array_ptr),
1030 defParamsArray(def_params_array_ptr), mField(m_field) {
1031 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR,
"-poisson_ratio",
1033 rval = mField.get_moab().tag_get_handle(
1034 "YOUNGS_MODULUS", 1, MB_TYPE_DOUBLE, thYoungModulus, MB_TAG_SPARSE);
1035 if (rval != MB_SUCCESS) {
1036 MOFEM_LOG(
"ADOLCPlasticityWold", Sev::warning)
1037 <<
"YOUNGS_MODULUS tag does not exist using default value";
1038 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR,
"-young_modulus",
1042 rval = mField.get_moab().tag_get_handle(
1043 "Yield_Strength_C", 1, MB_TYPE_DOUBLE, thCompressiveYieldStress,
1046 if (rval != MB_SUCCESS) {
1047 MOFEM_LOG(
"ADOLCPlasticityWold", Sev::warning)
1048 <<
"Yield_Strength_C tag does not exist using default value";
1049 yieldStrengthC = defParamsArray->at(SIGMA_YC);
1052 rval = mField.get_moab().tag_get_handle(
1053 "Yield_Strength_T", 1, MB_TYPE_DOUBLE, thTensionYieldStress,
1056 if (rval != MB_SUCCESS) {
1057 MOFEM_LOG(
"ADOLCPlasticityWold", Sev::warning)
1058 <<
"Yield_Strength_T tag does not exist using default value";
1059 yieldStrengthT = defParamsArray->at(SIGMA_YT);
1063 "Can not get data from block");
1066 MoFEMErrorCode doWork(
int side, EntityType type,
1067 EntitiesFieldData::EntData &data) {
1070 for (
auto &b : blockData) {
1072 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
1075 if (thYoungModulus != 0)
1076 CHKERR mField.get_moab().tag_get_data(thYoungModulus, &fe_ent, 1,
1079 if (thCompressiveYieldStress != 0)
1080 CHKERR mField.get_moab().tag_get_data(
1081 thCompressiveYieldStress, &fe_ent, 1, &yieldStrengthC);
1083 if (thTensionYieldStress != 0)
1084 CHKERR mField.get_moab().tag_get_data(
1085 thTensionYieldStress, &fe_ent, 1, &yieldStrengthT);
1087 std::copy(b.paramsArray.begin(), b.paramsArray.end(),
1088 paramsArrayPtr->begin());
1090 (*paramsArrayPtr)[
LAMBDA] =
LAMBDA(youngsModulus, poissonRatio);
1091 (*paramsArrayPtr)[
MU] =
MU(youngsModulus, poissonRatio);
1092 (*paramsArrayPtr)[SIGMA_YC] = yieldStrengthC;
1093 (*paramsArrayPtr)[SIGMA_YT] = yieldStrengthT;
1099 std::copy(defParamsArray->begin(), defParamsArray->end(),
1100 paramsArrayPtr->begin());
1106 boost::shared_ptr<ParamsArray> paramsArrayPtr;
1107 boost::shared_ptr<ParamsArray> defParamsArray;
1115 std::vector<BlockData> blockData;
1119 std::vector<const CubitMeshSets *> meshset_vec_ptr,
1123 for (
auto m : meshset_vec_ptr) {
1125 "Heterogeneous Paraboloidal")
1127 std::vector<double> block_data;
1128 CHKERR m->getAttributes(block_data);
1129 if (block_data.size() != 0) {
1131 "Expected that block has zero attributes");
1133 auto get_block_ents = [&]() {
1136 m_field.
get_moab().get_entities_by_handle(
m->meshset, ents,
true);
1140 blockData.push_back({*defParamsArray, get_block_ents()});
1141 std::copy(block_data.begin(), block_data.end(),
1142 blockData.back().paramsArray.begin());
1150 <<
"Young modulus: " << (blockData.back().paramsArray)[
LAMBDA];
1152 <<
"Poisson ratio: " << (blockData.back().paramsArray)[
MU];
1154 blockData.back().paramsArray[
LAMBDA] =
1159 <<
"LAMBDA: " << (blockData.back().paramsArray)[
LAMBDA];
1161 <<
"MU: " << (blockData.back().paramsArray)[
MU];
1163 <<
"NUP: " << (blockData.back().paramsArray)[NUP];
1165 <<
"SIGMA_YT: " << (blockData.back().paramsArray)[SIGMA_YT];
1167 <<
"SIGMA_YC: " << (blockData.back().paramsArray)[SIGMA_YC];
1169 <<
"HT: " << (blockData.back().paramsArray)[HT];
1171 <<
"HC: " << (blockData.back().paramsArray)[HC];
1173 <<
"NT: " << (blockData.back().paramsArray)[NT];
1175 <<
"NC: " << (blockData.back().paramsArray)[NC];
1181 auto cubit_meshsets_vec_ptr =
1182 m_field.
getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
1183 std::regex((boost::format(
"%s(.*)") % block_name).str()));
1186 m_field, sev, cubit_meshsets_vec_ptr));