v0.16.3
Loading...
Searching...
No Matches
MatElastic.cpp
Go to the documentation of this file.
1/**
2 * @file MatElastic.cpp
3 * @author your name (you@domain.com)
4 * @brief
5 * @version 0.1
6 * @date 2026-03-26
7 *
8 * @copyright Copyright (c) 2026
9 *
10 */
11
12#include <MoFEM.hpp>
13
14using namespace MoFEM;
15
16#include "MatOps.hpp"
17#include "MatElastic.hpp"
18
19namespace MatOps {
20
21template <int DIM> struct MatPiolaResponseImpl : public MatPiolaResponse {
22 using MatPiolaResponse::MatPiolaResponse;
24 createOp(boost::shared_ptr<PhysicalEquations> physical_ptr, bool eval_stress,
25 bool eval_tangent, bool update) override;
26};
27
28template <int DIM> struct MatElasticImpl : public MatElastic {
29 using MatElastic::MatElastic;
31 createOp(boost::shared_ptr<PhysicalEquations> physical_ptr, bool eval_stress,
32 bool eval_tangent, bool update) override;
33};
34
35template <int DIM> struct TopoMatObjectiveImpl : public PhysicalEquations {
38 createOp(boost::shared_ptr<PhysicalEquations> physical_ptr, bool eval_stress,
39 bool eval_tangent, bool update) override;
40};
41
42template <int DIM, typename DomainEleOp>
44
45 using OP = DomainEleOp;
46
47 OpEvalMatPiolaResponseImpl(boost::shared_ptr<PhysicalEquations> physical_ptr,
48 bool eval_stress, bool eval_tangent, bool update)
49 : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), physicalPtr(physical_ptr),
50 evalStress(eval_stress), evalTangent(eval_tangent),
51 updateState(update) {}
52
53 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
54
55protected:
56 boost::shared_ptr<PhysicalEquations> physicalPtr;
57 const bool evalStress;
58 const bool evalTangent;
59 const bool updateState;
60};
61
62template <int DIM, typename DomainEleOp>
65 EntData &data) {
67
68 int nb_integration_pts = OP::getGaussPts().size2();
69
70 auto get_tag = [&]() {
71 if (physicalPtr->tagVsRangePtr) {
72 for (const auto &tag_range_pair : *(physicalPtr->tagVsRangePtr)) {
73 if (tag_range_pair.second.find(DomainEleOp::getFEEntityHandle()) !=
74 tag_range_pair.second.end()) {
75 return tag_range_pair.first;
76 }
77 }
78 }
79#ifndef NDEBUG
80 if (MatOpsTagsRegistry::getTagName(physicalPtr->tAg).empty()) {
82 "ADOL-C tag not found " +
83 std::to_string(physicalPtr->tAg));
84 }
85#endif
86 return physicalPtr->tAg; // Default tag if no range matches
87 };
88
89 const int current_tag = get_tag();
90 auto *fe_ptr = const_cast<FEMethod *>(this->getFEMethod());
91 const auto ent = fe_ptr->getFEEntityHandle();
92
93 physicalPtr->matOpsDataPtr->getActiveDataPtr("F")->resize(DIM, DIM, false);
94 physicalPtr->matOpsDataPtr->getDependentDataPtr("P")->resize(DIM, DIM, false);
95 physicalPtr->matOpsDataPtr->getDependentDerivativesDataPtr("P_dF")->resize(
96 DIM * DIM, DIM * DIM, false);
97 auto mat_grad_ptr = physicalPtr->matOpsDataPtr->getCommonDataPtr("grad");
98 auto mat_P_ptr = physicalPtr->matOpsDataPtr->getCommonDataPtr("P");
99 auto mat_P_dF_ptr = physicalPtr->matOpsDataPtr->getCommonDataPtr("P_dF");
100
101#ifndef NDEBUG
102 if (!mat_grad_ptr || !mat_P_ptr || !mat_P_dF_ptr) {
103 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
104 "Missing common data for ADOL-C evaluation");
105 }
106 if (mat_grad_ptr->size2() != DIM * DIM) {
107 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
108 "Inconsistent size of gradient matrix for ADOL-C evaluation");
109 }
110 if (mat_grad_ptr->size1() != nb_integration_pts) {
111 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
112 "Inconsistent size of gradient matrix data for ADOL-C evaluation "
113 "%zu != %d",
114 mat_grad_ptr->size1(), nb_integration_pts);
115 }
116#endif
117
119 auto get_grad_at_pts =
120 MatrixSizeHelper<GetFTensor2FromMatType<DIM, DIM, -1, DL>, DL>::get(
121 *mat_grad_ptr, nb_integration_pts);
122 auto get_P_at_pts =
124 *mat_P_ptr, nb_integration_pts);
125 auto get_P_dF_at_pts =
126 MatrixSizeHelper<GetFTensor4FromMatType<DIM, DIM, DIM, DIM, -1, DL>,
127 DL>::size(*mat_P_dF_ptr, nb_integration_pts);
128
129 if (evalStress) {
130 FTENSOR_INDEX(DIM, i);
131 FTENSOR_INDEX(DIM, J);
132
133 auto t_grad_at_pts = get_grad_at_pts();
134 auto t_P_at_pts = get_P_at_pts();
135
136 auto next = [&]() {
137 ++t_grad_at_pts;
138 ++t_P_at_pts;
139 };
140
141 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
142 auto t_F = getFTensor2FromPtr<DIM, DIM>(
143 physicalPtr->matOpsDataPtr->getActiveDataPtr("F")->data().data());
144 t_F(i, J) = t_grad_at_pts(i, J);
145 CHKERR physicalPtr->setParams(fe_ptr, gg);
146 CHKERR physicalPtr->evaluateVariable(current_tag, ent, gg);
147 auto t_P = getFTensor2FromPtr<DIM, DIM>(
148 physicalPtr->matOpsDataPtr->getDependentDataPtr("P")->data().data());
149 t_P_at_pts(i, J) = t_P(i, J);
150 next();
151 }
152 }
153
154 if (evalTangent) {
155 FTENSOR_INDEX(DIM, i);
156 FTENSOR_INDEX(DIM, J);
157 FTENSOR_INDEX(DIM, k);
158 FTENSOR_INDEX(DIM, L);
159
160 auto t_grad_at_pts = get_grad_at_pts();
161 auto t_P_dF_at_pts = get_P_dF_at_pts();
162 auto next = [&]() {
163 ++t_grad_at_pts;
164 ++t_P_dF_at_pts;
165 };
166
167 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
168 auto t_F = getFTensor2FromPtr<DIM, DIM>(
169 physicalPtr->matOpsDataPtr->getActiveDataPtr("F")->data().data());
170 t_F(i, J) = t_grad_at_pts(i, J);
171 CHKERR physicalPtr->setParams(fe_ptr, gg);
172 CHKERR physicalPtr->evaluateDerivatives(current_tag, ent, gg);
173 auto t_P_dF = getFTensor4FromPtr<DIM, DIM, DIM, DIM>(
174 physicalPtr->matOpsDataPtr->getDependentDerivativesDataPtr("P_dF")
175 ->data()
176 .data());
177 t_P_dF_at_pts(i, J, k, L) = t_P_dF(i, J, k, L);
178 next();
179 }
180 }
181
182 if (updateState) {
183 FTENSOR_INDEX(DIM, i);
184 FTENSOR_INDEX(DIM, J);
185
186 auto t_grad_at_pts = get_grad_at_pts();
187 auto next = [&]() { ++t_grad_at_pts; };
188
189 for (int gg = 0; gg != nb_integration_pts; ++gg) {
190 auto t_F = getFTensor2FromPtr<DIM, DIM>(
191 physicalPtr->matOpsDataPtr->getActiveDataPtr("F")->data().data());
192 t_F(i, J) = t_grad_at_pts(i, J);
193 CHKERR physicalPtr->setParams(fe_ptr, gg);
194 CHKERR physicalPtr->updateState(current_tag, ent, gg);
195 next();
196 }
197 }
198
200}
201
202template <int DIM, typename DomainEleOp>
204
206
207 OpEvalTopoMatObjectiveImpl(boost::shared_ptr<PhysicalEquations> physical_ptr,
208 bool eval_stress, bool eval_tangent, bool update)
209 : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), physicalPtr(physical_ptr),
210 evalStress(eval_stress), evalTangent(eval_tangent),
211 updateState(update) {}
212
213 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
214
215protected:
216 boost::shared_ptr<PhysicalEquations> physicalPtr;
217 const bool evalStress;
218 const bool evalTangent;
219 const bool updateState;
220};
221
222template <int DIM, typename DomainEleOp>
225 EntData &data) {
227
228 int nb_integration_pts = OP::getGaussPts().size2();
229
230 auto get_tag = [&]() {
231 if (physicalPtr->tagVsRangePtr) {
232 for (const auto &tag_range_pair : *(physicalPtr->tagVsRangePtr)) {
233 if (tag_range_pair.second.find(DomainEleOp::getFEEntityHandle()) !=
234 tag_range_pair.second.end()) {
235 return tag_range_pair.first;
236 }
237 }
238 }
239#ifndef NDEBUG
240 if (MatOpsTagsRegistry::getTagName(physicalPtr->tAg).empty()) {
242 "ADOL-C tag not found " +
243 std::to_string(physicalPtr->tAg));
244 }
245#endif
246 return physicalPtr->tAg; // Default tag if no range matches
247 };
248
249 const int current_tag = get_tag();
250 auto *fe_ptr = const_cast<FEMethod *>(this->getFEMethod());
251 const auto ent = fe_ptr->getFEEntityHandle();
252
253 physicalPtr->matOpsDataPtr->getActiveDataPtr("Jac")->resize(DIM, DIM, false);
254 physicalPtr->matOpsDataPtr->getDependentDataPtr("Obj")->resize(1, false);
255 physicalPtr->matOpsDataPtr->getDependentDerivativesDataPtr("dObj")->resize(
256 DIM, DIM, false);
257 auto mat_jac_ptr = physicalPtr->matOpsDataPtr->getCommonDataPtr("Jac");
258 auto mat_grad_ptr = physicalPtr->matOpsDataPtr->getCommonDataPtr("grad");
259 auto mat_var_grad_ptr =
260 physicalPtr->matOpsDataPtr->getCommonDataPtr("var_grad");
261 auto mat_obj_ptr = physicalPtr->matOpsDataPtr->getCommonDataPtr("Obj");
262 auto mat_dObj_ptr = physicalPtr->matOpsDataPtr->getCommonDataPtr("dObj");
263
264#ifndef NDEBUG
265 if (!mat_jac_ptr || !mat_grad_ptr || !mat_var_grad_ptr || !mat_obj_ptr ||
266 !mat_dObj_ptr) {
267 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
268 "Missing common data for ADOL-C evaluation");
269 }
270#endif
271
273 auto get_jac_at_pts =
274 MatrixSizeHelper<GetFTensor2FromMatType<DIM, DIM, -1, DL>, DL>::get(
275 *mat_jac_ptr, nb_integration_pts);
276 [[maybe_unused]] auto get_grad_at_pts =
277 MatrixSizeHelper<GetFTensor2FromMatType<DIM, DIM, -1, DL>, DL>::get(
278 *mat_grad_ptr, nb_integration_pts);
279 [[maybe_unused]] auto get_var_grad_at_pts =
280 MatrixSizeHelper<GetFTensor2FromMatType<DIM, DIM, -1, DL>, DL>::get(
281 *mat_var_grad_ptr, nb_integration_pts);
282 auto get_obj_at_pts =
284 *mat_obj_ptr, nb_integration_pts);
285 auto get_dObj_at_pts =
287 *mat_dObj_ptr, nb_integration_pts);
288
289 if (evalStress) {
290 FTENSOR_INDEX(DIM, i);
291 FTENSOR_INDEX(DIM, J);
292
293 auto t_jac_at_pts = get_jac_at_pts();
294 auto t_obj_at_pts = get_obj_at_pts();
295
296 auto next = [&]() {
297 ++t_jac_at_pts;
298 ++t_obj_at_pts;
299 };
300
301 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
302 auto t_jac = getFTensor2FromPtr<DIM, DIM>(
303 physicalPtr->matOpsDataPtr->getActiveDataPtr("Jac")->data().data());
304 t_jac(i, J) = t_jac_at_pts(i, J);
305 CHKERR physicalPtr->setParams(fe_ptr, gg);
306 CHKERR physicalPtr->evaluateVariable(current_tag, ent, gg);
307 auto t_obj = getFTensor1FromPtr<1, 1>(
308 physicalPtr->matOpsDataPtr->getDependentDataPtr("Obj")
309 ->data()
310 .data());
311 t_obj_at_pts(0) = t_obj(0);
312 next();
313 }
314 }
315
316 if (evalTangent) {
317 FTENSOR_INDEX(DIM, i);
318 FTENSOR_INDEX(DIM, J);
319 FTENSOR_INDEX(DIM, k);
320 FTENSOR_INDEX(DIM, L);
321
322 auto t_jac_at_pts = get_jac_at_pts();
323 auto t_dObj_at_pts = get_dObj_at_pts();
324 auto next = [&]() {
325 ++t_jac_at_pts;
326 ++t_dObj_at_pts;
327 };
328
329 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
330 auto t_jac = getFTensor2FromPtr<DIM, DIM>(
331 physicalPtr->matOpsDataPtr->getActiveDataPtr("Jac")->data().data());
332 t_jac(i, J) = t_jac_at_pts(i, J);
333 CHKERR physicalPtr->setParams(fe_ptr, gg);
334 CHKERR physicalPtr->evaluateDerivatives(current_tag, ent, gg);
335 auto t_dObj = getFTensor2FromPtr<DIM, DIM>(
336 physicalPtr->matOpsDataPtr->getDependentDerivativesDataPtr("dObj")
337 ->data()
338 .data());
339 t_dObj_at_pts(i, J) = t_dObj(i, J);
340 next();
341 }
342 }
343
344 if (updateState) {
345 }
346
348}
349
350template <int DIM, typename EleOp>
351static EleOp *createOpImpl(boost::shared_ptr<PhysicalEquations> physical_ptr,
352 bool eval_stress, bool eval_tangent, bool update) {
353 return new OpEvalMatPiolaResponseImpl<DIM, EleOp>(physical_ptr, eval_stress,
354 eval_tangent, update);
355}
356
357template <>
359MatElasticImpl<3>::createOp(boost::shared_ptr<PhysicalEquations> physical_ptr,
360 bool eval_stress, bool eval_tangent, bool update) {
362 return createOpImpl<3, EleOp>(physical_ptr, eval_stress, eval_tangent,
363 update);
364}
365
366template <>
368 boost::shared_ptr<PhysicalEquations> physical_ptr, bool eval_stress,
369 bool eval_tangent, bool update) {
371 return createOpImpl<3, EleOp>(physical_ptr, eval_stress, eval_tangent,
372 update);
373}
374
375template <>
377 boost::shared_ptr<PhysicalEquations> physical_ptr, bool eval_stress,
378 bool eval_tangent, bool update) {
380 return createOpImpl<2, EleOp>(physical_ptr, eval_stress, eval_tangent,
381 update);
382}
383
384template <>
386MatElasticImpl<2>::createOp(boost::shared_ptr<PhysicalEquations> physical_ptr,
387 bool eval_stress, bool eval_tangent, bool update) {
389 return createOpImpl<2, EleOp>(physical_ptr, eval_stress, eval_tangent,
390 update);
391}
392
393template <int DIM, typename EleOp>
394static EleOp *
395createTopoOpImpl(boost::shared_ptr<PhysicalEquations> physical_ptr,
396 bool eval_stress, bool eval_tangent, bool update) {
397 return new OpEvalTopoMatObjectiveImpl<DIM, EleOp>(physical_ptr, eval_stress,
398 eval_tangent, update);
399}
400
401template <>
403 boost::shared_ptr<PhysicalEquations> physical_ptr, bool eval_stress,
404 bool eval_tangent, bool update) {
406 return createTopoOpImpl<3, EleOp>(physical_ptr, eval_stress, eval_tangent,
407 update);
408}
409
410template <>
412 boost::shared_ptr<PhysicalEquations> physical_ptr, bool eval_stress,
413 bool eval_tangent, bool update) {
415 return createTopoOpImpl<2, EleOp>(physical_ptr, eval_stress, eval_tangent,
416 update);
417}
418
419} // namespace MatOps
420
421#include "MatNeohookean.cpp"
426#include "MatGenericElastic.cpp"
std::string type
Dispatch first-Piola responses without implying stored energy.
Abaqus-compatible compressible Neo-Hookean material.
ADOL-C implementation of the volume-length quality material.
PetscInt size
#define FTENSOR_INDEX(DIM, I)
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
@ NOSPACE
Definition definitions.h:83
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'J', DIM1 > J
Definition level_set.cpp:30
FTensor::Index< 'k', 3 > k
static EleOp * createTopoOpImpl(boost::shared_ptr< PhysicalEquations > physical_ptr, bool eval_stress, bool eval_tangent, bool update)
static EleOp * createOpImpl(boost::shared_ptr< PhysicalEquations > physical_ptr, bool eval_stress, bool eval_tangent, bool update)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
decltype(GetFTensor4FromMatImpl< Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, Tensor_Dim3, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor4FromMatType
decltype(GetFTensor1FromMatImpl< Tensor_Dim, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor1FromMatType
decltype(GetFTensor2FromMatImpl< Tensor_Dim0, Tensor_Dim1, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor2FromMatType
ForcesAndSourcesCore::UserDataOperator * createOp(boost::shared_ptr< PhysicalEquations > physical_ptr, bool eval_stress, bool eval_tangent, bool update) override
Conservative elastic response with a stored-energy density.
static std::string getTagName(int tag)
Definition MatOps.cpp:90
ForcesAndSourcesCore::UserDataOperator * createOp(boost::shared_ptr< PhysicalEquations > physical_ptr, bool eval_stress, bool eval_tangent, bool update) override
Constitutive response expressed by first Piola stress.
Definition MatOps.hpp:186
OpEvalMatPiolaResponseImpl(boost::shared_ptr< PhysicalEquations > physical_ptr, bool eval_stress, bool eval_tangent, bool update)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< PhysicalEquations > physicalPtr
boost::shared_ptr< PhysicalEquations > physicalPtr
OpEvalTopoMatObjectiveImpl(boost::shared_ptr< PhysicalEquations > physical_ptr, bool eval_stress, bool eval_tangent, bool update)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
ForcesAndSourcesCore::UserDataOperator * createOp(boost::shared_ptr< PhysicalEquations > physical_ptr, bool eval_stress, bool eval_tangent, bool update) override
Data on single entity (This is passed as argument to DataOperator::doWork)
Structure for user loop methods on finite elements.
EntityHandle getFEEntityHandle() const
Get the entity handle of the current finite element.