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,
63 MoFEM::Interface &m_field, Sev sev,
64 std::vector<const CubitMeshSets *> meshset_vec_ptr,
65 double scale)
66 : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), matDPtr(m),
67 bulkModulusKDefault(bulk_modulus_K),
68 shearModulusGDefault(shear_modulus_G), isPlaneStrain(is_plane_strain),
69 scaleYoungModulus(scale) {
70 CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
71 "Can not get data from block");
72 }
73
74 MoFEMErrorCode doWork(int side, EntityType type,
75 EntitiesFieldData::EntData &data) {
77
78 for (auto &b : blockData) {
79
80 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
81 CHKERR getMatDPtr(matDPtr, b.bulkModulusK * scaleYoungModulus,
82 b.shearModulusG * scaleYoungModulus, isPlaneStrain);
84 }
85 }
86
87 CHKERR getMatDPtr(matDPtr, bulkModulusKDefault * scaleYoungModulus,
88 shearModulusGDefault * scaleYoungModulus,
89 isPlaneStrain);
91 }
92
93 private:
94 boost::shared_ptr<MatrixDouble> matDPtr;
95 const double scaleYoungModulus;
96
97 struct BlockData {
98 double bulkModulusK; ///< Bulk modulus
99 double shearModulusG; ///< Shear modulus
100 Range blockEnts; ///< Entities in the block
101 };
102
103 double bulkModulusKDefault;
104 double shearModulusGDefault;
105 std::vector<BlockData> blockData;
106 bool isPlaneStrain; ///< Plane strain flag
107
108 MoFEMErrorCode
109 extractBlockData(MoFEM::Interface &m_field,
110 std::vector<const CubitMeshSets *> meshset_vec_ptr,
111 Sev sev) {
113
114 for (auto m : meshset_vec_ptr) {
115 MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock") << *m;
116 std::vector<double> block_data;
117 CHKERR m->getAttributes(block_data);
118 if (block_data.size() < 2) {
119 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
120 "Expected that block has two attribute");
121 }
122 auto get_block_ents = [&]() {
123 Range ents;
124 CHKERR
125 m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
126 return ents;
127 };
128
129 double young_modulus = block_data[0];
130 double poisson_ratio = block_data[1];
131 double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
132 double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
133
134 MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
135 << "E = " << young_modulus << " nu = " << poisson_ratio;
136
137 blockData.push_back(
138 {bulk_modulus_K, shear_modulus_G, get_block_ents()});
139 }
140 MOFEM_LOG_CHANNEL("WORLD");
142 }
143 //! [Calculate elasticity tensor]
144 // Compute elasticity tensor D
145 MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
146 double bulk_modulus_K, double shear_modulus_G,
147 bool is_plane_strain) {
149
150 auto set_material_stiffness = [&]() {
151 FTensor::Index<'i', DIM> i;
152 FTensor::Index<'j', DIM> j;
153 FTensor::Index<'k', DIM> k;
154 FTensor::Index<'l', DIM> l;
156 double A = 1.;
157 /**
158 * @brief Toggle for plane strain.
159 *
160 * @note If `DIM == 2` and `is_plane_strain` is true, the code uses the
161 * plane strain elasticity tensor.
162 * @note If `DIM == 2` and `is_plane_strain` is false, the code uses the
163 * plane stress elasticity tensor.
164 */
165 if (DIM == 2 && !is_plane_strain) {
166 A = 2 * shear_modulus_G /
167 (bulk_modulus_K + (4. / 3.) * shear_modulus_G);
168 }
169 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mat_D_ptr);
170 t_D(i, j, k, l) =
171 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
172 A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) * t_kd(i, j) *
173 t_kd(k, l);
174 };
175
176 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
177 mat_D_ptr->resize(size_symm * size_symm, 1);
178 set_material_stiffness();
180 }
181 //! [Calculate elasticity tensor]
182 };
183
184 double E = 1.0;
185 double nu = 0.3;
186
187 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-young_modulus", &E,
188 PETSC_NULLPTR);
189 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-poisson_ratio", &nu,
190 PETSC_NULLPTR);
191
192 PetscBool is_plane_strain = PETSC_FALSE;
193 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-plane_strain",
194 &is_plane_strain, PETSC_NULLPTR);
195
196 double bulk_modulus_K = E / (3 * (1 - 2 * nu));
197 double shear_modulus_G = E / (2 * (1 + nu));
198 pip.push_back(new OpMatBlocks(
199 mat_D_Ptr, bulk_modulus_K, shear_modulus_G, is_plane_strain, m_field, sev,
200
201 // Get blockset using regular expression
202 m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
203
204 (boost::format("%s(.*)") % block_name).str()
205
206 )),
207 scale
208
209 ));
210
212}
213// Factory to initialize common data to get strain, stress and D
214template <int DIM, IntegrationType I, typename DomainEleOp>
216 MoFEM::Interface &m_field,
217 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
218 std::string field_name, std::string block_name, Sev sev, double scale = 1) {
219
220 auto common_ptr = boost::make_shared<HookeOps::CommonData>();
221 common_ptr->matDPtr = boost::make_shared<MatrixDouble>();
222 common_ptr->matGradPtr = boost::make_shared<MatrixDouble>();
223
224 CHK_THROW_MESSAGE(addMatBlockOps<DIM>(m_field, pip, block_name,
225 common_ptr->matDPtr, sev, scale),
226 "addMatBlockOps");
227
228 pip.push_back(new OpCalculateVectorFieldGradient<DIM, DIM>(
229 field_name, common_ptr->matGradPtr));
230 pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(common_ptr->matGradPtr,
231 common_ptr->getMatStrain()));
232 pip.push_back(new OpTensorTimesSymmetricTensor<SPACE_DIM, SPACE_DIM>(
233 common_ptr->getMatStrain(), common_ptr->getMatCauchyStress(),
234 common_ptr->matDPtr));
235
236 return common_ptr;
237}
238
239/**
240@brief Factory function to create and push internal force operators for the
241domain RHS. (RHS first overload)
242 * Assemble domain right-hand-side vector- internal force vector.
243 * @param m_field MoFEM interface reference
244 * @param field_name Field name string
245 * @param common_ptr Hookeian Material data pointer
246 * @param is_non_linear option for solving linear nonlinear
247 */
248template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
249MoFEMErrorCode opFactoryDomainRhs(
250 MoFEM::Interface &m_field,
251 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
252 std::string field_name, boost::shared_ptr<HookeOps::CommonData> common_ptr,
253 Sev sev, bool is_non_linear = false) {
255 //! [Push Internal forces]
256 using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
257 A>::template LinearForm<I>;
258 using OpInternalForce =
259 typename B::template OpGradTimesSymTensor<1, DIM, DIM>;
260
261 pip.push_back(
262 new OpInternalForce(field_name, common_ptr->getMatCauchyStress(),
263 [is_non_linear](double, double, double) constexpr {
264 return (is_non_linear ? 1 : -1);
265 }));
266 //! [Push Internal forces]
268}
269
270/**
271 * @brief Overloaded factory for domain RHS for material scaling (RHS Second overload)
272 *
273 * Initializes and pushes operators to assemble the domain RHS vector,
274 * representing internal forces, with optional scaling.
275 * @tparam DIM Spatial dimension
276 * @tparam A Assembly type
277 * @tparam I Integration type
278
279 * @param pip Operator container
280 * @param field_name Field variable string
281 * @param block_name Mesh block name
282 * @param sev Logging severity level
283 * @param is_non_linear Linear or Nonlinear problem flag
284 * @param scale Yonuns modulas scaling factor
285
286 */
287template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
288MoFEMErrorCode opFactoryDomainRhs(
289 MoFEM::Interface &m_field,
290 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
291 std::string field_name, std::string block_name, Sev sev,
292 bool is_non_linear = false, double scale = 1) {
294
295 auto common_ptr = commonDataFactory<DIM, I, DomainEleOp>(
296 m_field, pip, field_name, block_name, sev, scale);
297 CHKERR opFactoryDomainRhs<DIM, A, I, DomainEleOp>(
298 m_field, pip, field_name, common_ptr, sev, is_non_linear);
299
301}
302/**
303 * @brief Assemble domain LHS K factory (LHS first overload)
304 * Initializes and pushes operators to assemble the LHS
305
306 */
307
308template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
309MoFEMErrorCode opFactoryDomainLhs(
310 MoFEM::Interface &m_field,
311 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
312 std::string field_name, boost::shared_ptr<HookeOps::CommonData> common_ptr,
313 Sev sev) {
315 //! [OpK]
316 using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
317 A>::template BiLinearForm<I>;
318
319 using OpKCauchy = typename B::template OpGradSymTensorGrad<1, DIM, DIM, 0>;
320 //! [OpK]
321 // Assumes constant D matrix per entity
322 pip.push_back(new OpKCauchy(field_name, field_name, common_ptr->matDPtr));
323
325}
326/**
327 * @brief Overloaded LHS factory with optional scaling (LHS Second overload)
328 *
329 * @param scale Young’s modulus scale
330 */
331
332template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
333MoFEMErrorCode opFactoryDomainLhs(
334 MoFEM::Interface &m_field,
335 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
336 std::string field_name, std::string block_name, Sev sev, double scale = 1) {
338
339 auto common_ptr = commonDataFactory<DIM, I, DomainEleOp>(
340 m_field, pip, field_name, block_name, sev, scale);
341 CHKERR opFactoryDomainLhs<DIM, A, I, DomainEleOp>(m_field, pip, field_name,
342 common_ptr, sev);
343
345}
346} // namespace HookeOps
347
348#endif // __HOOKE_OPS_HPP__
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
PetscBool is_plane_strain
Definition adjoint.cpp:42
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()
@ 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 shear_modulus_G
double bulk_modulus_K
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:309
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:215
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:249
constexpr AssemblyType A
double scale
Definition plastic.cpp:123
double poisson_ratio
Poisson ratio.
Definition plastic.cpp:126
constexpr auto size_symm
Definition plastic.cpp:42
double young_modulus
Young modulus.
Definition plastic.cpp:125
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:49
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.