v0.15.0
Loading...
Searching...
No Matches
HookeOps.hpp
Go to the documentation of this file.
1/**
2 * @file HookeOps.hpp
3 * @brief Implementation of Hookes operator Hookes for linear elastic problems
4 * in MoFEM.
5 * @version 0.14.0
6 * @date 2025-04-10
7 */
8
9#ifndef __HOOKE_OPS_HPP__
10#define __HOOKE_OPS_HPP__
11
12namespace HookeOps {
13// Shared data structure for sotting .
14struct CommonData : public boost::enable_shared_from_this<CommonData> {
15 boost::shared_ptr<MatrixDouble> matGradPtr;
16 boost::shared_ptr<MatrixDouble> matDPtr;
17 MatrixDouble matCauchyStressPtr;
18 MatrixDouble matStrainPtr;
19
20 /**
21 * @brief Get shared Cauchy stress
22 * @return Shared Cauchy stress matrix
23 */
24 inline auto getMatCauchyStress() {
25 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
27 }
28 /**
29 * @brief Get shared strain
30 * @return Shared strain matrix
31 */
32 inline auto getMatStrain() {
33 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matStrainPtr);
34 }
35};
36
37/**
38 * @brief // Add operator to compute material stiffness tensor D based on block
39 * attributes.
40 *
41 * Extracts block-wise material properties and initializes elasticity tensor
42 * for each finite element block.
43 *
44 * @tparam DIM Problem dimension
45 * @param pip Operator pipeline container
46 * @param m_field MoFEM interface reference
47 * @param mat_D_Ptr Elasticity tensor pointer
48 * @param sev Logging Verbosity level
49 * @param scale Young's modulus scale
50 */
51template <int DIM>
52MoFEMErrorCode
54 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
55 std::string block_name,
56 boost::shared_ptr<MatrixDouble> mat_D_Ptr, Sev sev,
57 double scale = 1) {
59
60 struct OpMatBlocks : public DomainEleOp {
61 OpMatBlocks(boost::shared_ptr<MatrixDouble> m, double bulk_modulus_K,
62 double shear_modulus_G, MoFEM::Interface &m_field, Sev sev,
63 std::vector<const CubitMeshSets *> meshset_vec_ptr,
64 double scale)
65 : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), matDPtr(m),
66 bulkModulusKDefault(bulk_modulus_K),
67 shearModulusGDefault(shear_modulus_G), scaleYoungModulus(scale) {
68 CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
69 "Can not get data from block");
70 }
71
72 MoFEMErrorCode doWork(int side, EntityType type,
73 EntitiesFieldData::EntData &data) {
75
76 for (auto &b : blockData) {
77
78 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
79 CHKERR getMatDPtr(matDPtr, b.bulkModulusK * scaleYoungModulus,
80 b.shearModulusG * scaleYoungModulus);
82 }
83 }
84
85 CHKERR getMatDPtr(matDPtr, bulkModulusKDefault * scaleYoungModulus,
86 shearModulusGDefault * scaleYoungModulus);
88 }
89
90 private:
91 boost::shared_ptr<MatrixDouble> matDPtr;
92 const double scaleYoungModulus;
93
94 struct BlockData {
95 double bulkModulusK; ///< Bulk modulus
96 double shearModulusG; ///< Shear modulus
97 Range blockEnts; ///< Entities in the block
98 };
99
100 double bulkModulusKDefault;
101 double shearModulusGDefault;
102 std::vector<BlockData> blockData;
103
104 MoFEMErrorCode
105 extractBlockData(MoFEM::Interface &m_field,
106 std::vector<const CubitMeshSets *> meshset_vec_ptr,
107 Sev sev) {
109
110 for (auto m : meshset_vec_ptr) {
111 MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock") << *m;
112 std::vector<double> block_data;
113 CHKERR m->getAttributes(block_data);
114 if (block_data.size() < 2) {
115 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
116 "Expected that block has two attribute");
117 }
118 auto get_block_ents = [&]() {
119 Range ents;
120 CHKERR
121 m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
122 return ents;
123 };
124
125 double young_modulus = block_data[0];
126 double poisson_ratio = block_data[1];
127 double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
128 double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
129
130 MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
131 << "E = " << young_modulus << " nu = " << poisson_ratio;
132
133 blockData.push_back(
134 {bulk_modulus_K, shear_modulus_G, get_block_ents()});
135 }
136 MOFEM_LOG_CHANNEL("WORLD");
138 }
139 //! [Calculate elasticity tensor]
140 // Compute elasticity tensor D
141 MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
142 double bulk_modulus_K, double shear_modulus_G) {
144
145 auto set_material_stiffness = [&]() {
146 FTensor::Index<'i', DIM> i;
147 FTensor::Index<'j', DIM> j;
148 FTensor::Index<'k', DIM> k;
149 FTensor::Index<'l', DIM> l;
151 double A = 1.;
152 /**
153 * @brief Toggle for plane strain.
154 *
155 * @note If `DIM == 2` and `is_plane_strain` is true, the code uses the
156 * plane strain elasticity tensor.
157 * @note If `DIM == 2` and `is_plane_strain` is false, the code uses the
158 * plane stress elasticity tensor.
159 */
160 if (DIM == 2 && !is_plane_strain) {
161 A = 2 * shear_modulus_G /
162 (bulk_modulus_K + (4. / 3.) * shear_modulus_G);
163 }
164 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mat_D_ptr);
165 t_D(i, j, k, l) =
166 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
167 A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) * t_kd(i, j) *
168 t_kd(k, l);
169 };
170
171 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
172 mat_D_ptr->resize(size_symm * size_symm, 1);
173 set_material_stiffness();
175 }
176 //! [Calculate elasticity tensor]
177 };
178
179 double E = 1.0;
180 double nu = 0.3;
181
182 PetscOptionsBegin(PETSC_COMM_WORLD, "", "", "none");
183 CHKERR PetscOptionsScalar("-young_modulus", "Young modulus", "", E, &E,
184 PETSC_NULLPTR);
185 CHKERR PetscOptionsScalar("-poisson_ratio", "poisson ratio", "", nu, &nu,
186 PETSC_NULLPTR);
187 PetscOptionsEnd();
188
189 double bulk_modulus_K = E / (3 * (1 - 2 * nu));
190 double shear_modulus_G = E / (2 * (1 + nu));
191 pip.push_back(new OpMatBlocks(
192 mat_D_Ptr, bulk_modulus_K, shear_modulus_G, m_field, sev,
193
194 // Get blockset using regular expression
195 m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
196
197 (boost::format("%s(.*)") % block_name).str()
198
199 )),
200 scale
201
202 ));
203
205}
206// Factory to initialize common data to get strain, stress and D
207template <int DIM, IntegrationType I, typename DomainEleOp>
209 MoFEM::Interface &m_field,
210 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
211 std::string field_name, std::string block_name, Sev sev, double scale = 1) {
212
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>();
216
217 CHK_THROW_MESSAGE(addMatBlockOps<DIM>(m_field, pip, block_name,
218 common_ptr->matDPtr, sev, scale),
219 "addMatBlockOps");
220
221 pip.push_back(new OpCalculateVectorFieldGradient<DIM, DIM>(
222 field_name, common_ptr->matGradPtr));
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));
228
229 return common_ptr;
230}
231
232/**
233@brief Factory function to create and push internal force operators for the
234domain RHS. (RHS first overload)
235 * Assemble domain right-hand-side vector- internal force vector.
236 * @param m_field MoFEM interface reference
237 * @param field_name Field name string
238 * @param common_ptr Hookeian Material data pointer
239 * @param is_non_linear option for solving linear nonlinear
240 */
241template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
242MoFEMErrorCode opFactoryDomainRhs(
243 MoFEM::Interface &m_field,
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) {
248 //! [Push Internal forces]
249 using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
250 A>::template LinearForm<I>;
251 using OpInternalForce =
252 typename B::template OpGradTimesSymTensor<1, DIM, DIM>;
253
254 pip.push_back(
255 new OpInternalForce(field_name, common_ptr->getMatCauchyStress(),
256 [is_non_linear](double, double, double) constexpr {
257 return (is_non_linear ? 1 : -1);
258 }));
259 //! [Push Internal forces]
261}
262
263/**
264 * @brief Overloaded factory for domain RHS for material scaling (RHS Second overload)
265 *
266 * Initializes and pushes operators to assemble the domain RHS vector,
267 * representing internal forces, with optional scaling.
268 * @tparam DIM Spatial dimension
269 * @tparam A Assembly type
270 * @tparam I Integration type
271
272 * @param pip Operator container
273 * @param field_name Field variable string
274 * @param block_name Mesh block name
275 * @param sev Logging severity level
276 * @param is_non_linear Linear or Nonlinear problem flag
277 * @param scale Yonuns modulas scaling factor
278
279 */
280template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
281MoFEMErrorCode opFactoryDomainRhs(
282 MoFEM::Interface &m_field,
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) {
287
289 m_field, pip, field_name, block_name, sev, scale);
291 m_field, pip, field_name, common_ptr, sev, is_non_linear);
292
294}
295/**
296 * @brief Assemble domain LHS K factory (LHS first overload)
297 * Initializes and pushes operators to assemble the LHS
298
299 */
300
301template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
302MoFEMErrorCode opFactoryDomainLhs(
303 MoFEM::Interface &m_field,
304 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
305 std::string field_name, boost::shared_ptr<HookeOps::CommonData> common_ptr,
306 Sev sev) {
308 //! [OpK]
309 using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
310 A>::template BiLinearForm<I>;
311
312 using OpKCauchy = typename B::template OpGradSymTensorGrad<1, DIM, DIM, 0>;
313 //! [OpK]
314 // Assumes constant D matrix per entity
315 pip.push_back(new OpKCauchy(field_name, field_name, common_ptr->matDPtr));
316
318}
319/**
320 * @brief Overloaded LHS factory with optional scaling (LHS Second overload)
321 *
322 * @param scale Young’s modulus scale
323 */
324
325template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
326MoFEMErrorCode opFactoryDomainLhs(
327 MoFEM::Interface &m_field,
328 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
329 std::string field_name, std::string block_name, Sev sev, double scale = 1) {
331
333 m_field, pip, field_name, block_name, sev, scale);
335 common_ptr, sev);
336
338}
339} // namespace HookeOps
340
341#endif // __HOOKE_OPS_HPP__
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
Kronecker Delta class symmetric.
PetscBool is_plane_strain
Definition contact.cpp:95
#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()
@ NOSPACE
Definition definitions.h:83
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
double bulk_modulus_K
double shear_modulus_G
constexpr auto t_kd
#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.
Definition HookeOps.hpp:53
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...
Definition HookeOps.hpp:302
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)
Definition HookeOps.hpp:208
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)...
Definition HookeOps.hpp:242
double young_modulus
Young modulus.
Definition plastic.cpp:125
double poisson_ratio
Poisson ratio.
Definition plastic.cpp:126
double scale
Definition plastic.cpp:123
constexpr auto size_symm
Definition plastic.cpp:42
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)]
Definition seepage.cpp:48
FTensor::Index< 'm', 3 > m
boost::shared_ptr< MatrixDouble > matDPtr
Definition HookeOps.hpp:16
auto getMatStrain()
Get shared strain.
Definition HookeOps.hpp:32
auto getMatCauchyStress()
Get shared Cauchy stress.
Definition HookeOps.hpp:24
MatrixDouble matStrainPtr
Definition HookeOps.hpp:18
boost::shared_ptr< MatrixDouble > matGradPtr
Definition HookeOps.hpp:15
MatrixDouble matCauchyStressPtr
Definition HookeOps.hpp:17
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.