9#ifndef __HOOKE_OPS_HPP__
10#define __HOOKE_OPS_HPP__
14struct CommonData :
public boost::enable_shared_from_this<CommonData> {
25 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
33 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
matStrainPtr);
54 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
55 std::string block_name,
56 boost::shared_ptr<MatrixDouble> mat_D_Ptr, Sev sev,
63 std::vector<const CubitMeshSets *> meshset_vec_ptr,
69 "Can not get data from block");
72 MoFEMErrorCode doWork(
int side, EntityType type,
73 EntitiesFieldData::EntData &data) {
76 for (
auto &b : blockData) {
78 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
79 CHKERR getMatDPtr(matDPtr, b.bulkModulusK * scaleYoungModulus,
80 b.shearModulusG * scaleYoungModulus);
85 CHKERR getMatDPtr(matDPtr, bulkModulusKDefault * scaleYoungModulus,
86 shearModulusGDefault * scaleYoungModulus);
91 boost::shared_ptr<MatrixDouble> matDPtr;
92 const double scaleYoungModulus;
100 double bulkModulusKDefault;
101 double shearModulusGDefault;
102 std::vector<BlockData> blockData;
106 std::vector<const CubitMeshSets *> meshset_vec_ptr,
110 for (
auto m : meshset_vec_ptr) {
112 std::vector<double> block_data;
113 CHKERR m->getAttributes(block_data);
114 if (block_data.size() < 2) {
116 "Expected that block has two attribute");
118 auto get_block_ents = [&]() {
121 m_field.
get_moab().get_entities_by_handle(
m->meshset, ents,
true);
141 MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
145 auto set_material_stiffness = [&]() {
164 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mat_D_ptr);
171 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
173 set_material_stiffness();
182 PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"",
"none");
183 CHKERR PetscOptionsScalar(
"-young_modulus",
"Young modulus",
"",
E, &
E,
185 CHKERR PetscOptionsScalar(
"-poisson_ratio",
"poisson ratio",
"", nu, &nu,
191 pip.push_back(
new OpMatBlocks(
195 m_field.
getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
197 (boost::format(
"%s(.*)") % block_name).str()
207template <
int DIM, IntegrationType I,
typename DomainEleOp>
210 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
211 std::string
field_name, std::string block_name, Sev sev,
double scale = 1) {
213 auto common_ptr = boost::make_shared<HookeOps::CommonData>();
214 common_ptr->matDPtr = boost::make_shared<MatrixDouble>();
215 common_ptr->matGradPtr = boost::make_shared<MatrixDouble>();
218 common_ptr->matDPtr, sev,
scale),
221 pip.push_back(
new OpCalculateVectorFieldGradient<DIM, DIM>(
223 pip.push_back(
new OpSymmetrizeTensor<SPACE_DIM>(common_ptr->matGradPtr,
224 common_ptr->getMatStrain()));
225 pip.push_back(
new OpTensorTimesSymmetricTensor<SPACE_DIM, SPACE_DIM>(
226 common_ptr->getMatStrain(), common_ptr->getMatCauchyStress(),
227 common_ptr->matDPtr));
241template <
int DIM, AssemblyType A, IntegrationType I,
typename DomainEleOp>
244 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
245 std::string
field_name, boost::shared_ptr<HookeOps::CommonData> common_ptr,
246 Sev sev,
bool is_non_linear =
false) {
249 using B =
typename FormsIntegrators<DomainEleOp>::template Assembly<
251 using OpInternalForce =
252 typename B::template OpGradTimesSymTensor<1, DIM, DIM>;
255 new OpInternalForce(
field_name, common_ptr->getMatCauchyStress(),
256 [is_non_linear](
double,
double,
double)
constexpr {
257 return (is_non_linear ? 1 : -1);
280template <
int DIM, AssemblyType A, IntegrationType I,
typename DomainEleOp>
283 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
284 std::string
field_name, std::string block_name, Sev sev,
285 bool is_non_linear =
false,
double scale = 1) {
291 m_field, pip,
field_name, common_ptr, sev, is_non_linear);
301template <
int DIM, AssemblyType A, IntegrationType I,
typename DomainEleOp>
304 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
305 std::string
field_name, boost::shared_ptr<HookeOps::CommonData> common_ptr,
309 using B =
typename FormsIntegrators<DomainEleOp>::template Assembly<
312 using OpKCauchy =
typename B::template OpGradSymTensorGrad<1, DIM, DIM, 0>;
325template <
int DIM, AssemblyType A, IntegrationType I,
typename DomainEleOp>
328 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
329 std::string
field_name, std::string block_name, Sev sev,
double scale = 1) {
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
Kronecker Delta class symmetric.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
MoFEMErrorCode addMatBlockOps(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string block_name, boost::shared_ptr< MatrixDouble > mat_D_Ptr, Sev sev, double scale=1)
// Add operator to compute material stiffness tensor D based on block attributes.
MoFEMErrorCode opFactoryDomainLhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HookeOps::CommonData > common_ptr, Sev sev)
Assemble domain LHS K factory (LHS first overload) Initializes and pushes operators to assemble the L...
auto commonDataFactory(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, std::string block_name, Sev sev, double scale=1)
MoFEMErrorCode opFactoryDomainRhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HookeOps::CommonData > common_ptr, Sev sev, bool is_non_linear=false)
Factory function to create and push internal force operators for the domain RHS. (RHS first overload)...
double young_modulus
Young modulus.
double poisson_ratio
Poisson ratio.
constexpr auto field_name
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Only used with Hooke equation (linear material model)]
FTensor::Index< 'm', 3 > m
boost::shared_ptr< MatrixDouble > matDPtr
auto getMatStrain()
Get shared strain.
auto getMatCauchyStress()
Get shared Cauchy stress.
MatrixDouble matStrainPtr
boost::shared_ptr< MatrixDouble > matGradPtr
MatrixDouble matCauchyStressPtr
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.