v0.16.0
Loading...
Searching...
No Matches
MFrontInterface.cpp
Go to the documentation of this file.
1/** \file MFrontInterface.cpp
2 * \brief MFrontInterface
3 *
4 * MFrontInterface
5 *
6 */
7
8#include <MoFEM.hpp>
9#include <MFrontInterface.hpp>
10
11#ifdef WITH_MGIS
12
13 #include <MGIS/Behaviour/Behaviour.hxx>
14 #include <MGIS/Behaviour/BehaviourData.hxx>
15 #include "MGIS/LibrariesManager.hxx"
16 #include "MGIS/Behaviour/Integrate.hxx"
17
18using namespace FTensor;
19using namespace mgis::behaviour;
20
21namespace MoFEM {
22
23 constexpr double inv_sqr2 = boost::math::constants::half_root_two<double>();
24
25 #define VOIGT_VEC_SYMM_3D(VEC) \
26 VEC[0], inv_sqr2 *VEC[3], inv_sqr2 *VEC[4], VEC[1], inv_sqr2 *VEC[5], VEC[2]
27
28 #define VOIGT_VEC_SYMM_2D(VEC) VEC[0], inv_sqr2 *VEC[3], VEC[1]
29
30 #define VOIGT_VEC_SYMM_2D_FULL(VEC) \
31 VEC[0], inv_sqr2 *VEC[3], 0., VEC[1], 0., VEC[2]
32
33 #define VOIGT_VEC_3D(VEC) \
34 VEC[0], VEC[3], VEC[5], VEC[4], VEC[1], VEC[7], VEC[6], VEC[8], VEC[2]
35
36 #define VOIGT_VEC_2D(VEC) VEC[0], VEC[3], VEC[4], VEC[1]
37
38 #define VOIGT_VEC_2D_FULL(VEC) \
39 VEC[0], VEC[3], 0., VEC[4], VEC[1], 0., 0., 0., VEC[2]
40
41template <ModelHypothesis MH> struct MFrontEleType;
42
43template <> struct MFrontEleType<TRIDIMENSIONAL> {
44
45 MFrontEleType() = delete;
46 ~MFrontEleType() = delete;
47
48 using DomainEle = MoFEM::VolumeElementForcesAndSourcesCore;
49 using DomainEleOp = DomainEle::UserDataOperator;
50 using PostProcDomainEle = PostProcBrokenMeshInMoabBase<DomainEle>;
51
52 static constexpr int SPACE_DIM = 3;
53};
54
55template <> struct MFrontEleType<PLANESTRAIN> {
56
57 MFrontEleType() = delete;
58 ~MFrontEleType() = delete;
59
60 using DomainEle = MoFEM::FaceElementForcesAndSourcesCore;
61 using DomainEleOp = DomainEle::UserDataOperator;
62 using PostProcDomainEle = PostProcBrokenMeshInMoabBase<DomainEle>;
63
64 static constexpr int SPACE_DIM = 2;
65};
66
67template <> struct MFrontEleType<AXISYMMETRICAL> {
68
69 MFrontEleType() = delete;
70 ~MFrontEleType() = delete;
71
72 using DomainEle = FaceElementForcesAndSourcesCore;
73 using DomainEleOp = DomainEle::UserDataOperator;
74 using PostProcDomainEle = PostProcBrokenMeshInMoabBase<DomainEle>;
75
76 static constexpr int SPACE_DIM = 2;
77};
78
79
80template <ModelHypothesis MH, AssemblyType AT = AssemblyType::PETSC>
81struct MFrontInterfaceImpl : public MFrontInterface {
82
83 MFrontInterfaceImpl(MoFEM::Interface &m_field);
84
85 using DomainEle = typename MFrontEleType<MH>::DomainEle;
86 using DomainEleOp = typename MFrontEleType<MH>::DomainEleOp;
87 using PostProcDomainEle = typename MFrontEleType<MH>::PostProcDomainEle;
88
89 static constexpr int DIM = MFrontEleType<MH>::SPACE_DIM;
90
91 using OpInternalForce =
92 typename FormsIntegrators<DomainEleOp>::template Assembly<AT>::
93 template LinearForm<GAUSS>::template OpGradTimesTensor<1, DIM, DIM>;
94 using OpAssembleLhsFiniteStrains =
95 typename FormsIntegrators<DomainEleOp>::template Assembly<
96 AT>::template BiLinearForm<GAUSS>::template OpGradTensorGrad<1, DIM,
97 DIM, -1>;
98 using OpAssembleLhsSmallStrains =
99 typename FormsIntegrators<DomainEleOp>::template Assembly<AT>::
100 template BiLinearForm<GAUSS>::template OpGradSymTensorGrad<1, DIM,
101 DIM, 0>;
102
103 MoFEMErrorCode getCommandLineParameters() override;
104
105 MoFEMErrorCode opFactoryDomainRhs(
106 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
107 std::string field_name) override;
108
109 MoFEMErrorCode opFactoryDomainLhs(
110 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
111 std::string field_name) override;
112
113 MoFEMErrorCode
114 setUpdateInternalVariablesOperators(ForcesAndSourcesCore::RuleHookFun rule,
115 std::string field_name) override;
116
117 MoFEMErrorCode setPostProcessOperators(ForcesAndSourcesCore::RuleHookFun rule,
118 std::string fe_name,
119 std::string field_name,
120 int order) override;
121
122 MoFEMErrorCode updateInternalVariables(SmartPetscObj<DM> dm,
123 std::string fe_name) override;
124
125 MoFEMErrorCode postProcess(int step, SmartPetscObj<DM> dm,
126 string fe_name) override;
127
128 MoFEMErrorCode
129 setMonitor(boost::shared_ptr<MoFEM::FEMethod> monitor_ptr) override {
130 MoFEMFunctionBeginHot;
131 monitorPtr = monitor_ptr;
132 MoFEMFunctionReturnHot(0);
133 }
134
135private:
136 MoFEM::Interface &mField;
137
138 string optionsPrefix;
139
140 SmartPetscObj<DM> dM;
141
142 PetscBool saveGauss;
143 PetscBool saveDomain;
144
145 PetscBool testJacobian;
146 PetscReal randomFieldScale;
147
148 bool isFiniteKinematics;
149
150 FieldApproximationBase fieldBase;
151
152 boost::shared_ptr<PostProcDomainEle> postProcFe;
153 boost::shared_ptr<DomainEle> updateIntVariablesElePtr;
154
155 boost::shared_ptr<moab::Interface> moabGaussIntPtr;
156
157 boost::shared_ptr<MoFEM::FEMethod> monitorPtr;
158
159 boost::shared_ptr<CommonData> commonDataPtr;
160};
161
162template <ModelHypothesis MH>
163boost::shared_ptr<MFrontInterface>
164createMFrontInterfaceImpl(MoFEM::Interface &m_field, AssemblyType at) {
165
166 switch (at) {
167 case PETSC:
168 return boost::make_shared<MFrontInterfaceImpl<MH, PETSC>>(m_field);
169 case SCHUR:
170 return boost::make_shared<MFrontInterfaceImpl<MH, SCHUR>>(m_field);
171 case BLOCK_MAT:
172 return boost::make_shared<MFrontInterfaceImpl<MH, BLOCK_MAT>>(m_field);
173 case BLOCK_SCHUR:
174 return boost::make_shared<MFrontInterfaceImpl<MH, BLOCK_SCHUR>>(m_field);
175 default:
176 CHK_THROW_MESSAGE(MOFEM_NOT_IMPLEMENTED, "Unsupported AssemblyType");
177 }
178
179 return boost::shared_ptr<MFrontInterface>();
180}
181
182boost::shared_ptr<MFrontInterface>
183createMFrontInterface(MoFEM::Interface &m_field, ModelHypothesis mh,
184 AssemblyType at) {
185
186 switch (mh) {
187 case TRIDIMENSIONAL:
188 return createMFrontInterfaceImpl<TRIDIMENSIONAL>(m_field, at);
189 case PLANESTRAIN:
190 return createMFrontInterfaceImpl<PLANESTRAIN>(m_field, at);
191 case AXISYMMETRICAL:
192 return createMFrontInterfaceImpl<AXISYMMETRICAL>(m_field, at);
193 default:
194 CHK_THROW_MESSAGE(MOFEM_NOT_IMPLEMENTED, "Unsupported ModelHypothesis");
195 }
196 return boost::shared_ptr<MFrontInterface>();
197}
198
199struct MFrontInterface::CommonData {
200 CommonData(MoFEM::Interface &m_field) : mField(m_field) {};
201
202 MoFEMErrorCode setBlocks(int dim);
203 MoFEMErrorCode getInternalVar(const EntityHandle fe_ent,
204 const int nb_gauss_pts, const int var_size,
205 const int grad_size, const int stress_size,
206 bool is_large_strain = false);
207
208 MoFEMErrorCode setInternalVar(const EntityHandle fe_ent);
209 MoFEMErrorCode createTags();
210 MoFEMErrorCode clearTags();
211
212 MoFEM::Interface &mField;
213 boost::shared_ptr<MatrixDouble> mGradPtr;
214 boost::shared_ptr<MatrixDouble> mStressPtr;
215 boost::shared_ptr<MatrixDouble> mFullStrainPtr;
216 boost::shared_ptr<MatrixDouble> mFullStressPtr;
217
218 boost::shared_ptr<MatrixDouble> mPrevGradPtr;
219 boost::shared_ptr<MatrixDouble> mPrevStressPtr;
220
221 boost::shared_ptr<MatrixDouble> mDispPtr;
222 boost::shared_ptr<MatrixDouble> materialTangentPtr;
223 boost::shared_ptr<MatrixDouble> mFullTangentPtr;
224 boost::shared_ptr<MatrixDouble> internalVariablePtr;
225
226 struct BlockData;
227
228 std::map<int, BlockData> setOfBlocksData;
229 std::map<EntityHandle, int> blocksIDmap;
230
231 MatrixDouble lsMat;
232 MatrixDouble gaussPts;
233 MatrixDouble myN;
234
235 Tag internalVariableTag;
236 Tag stressTag;
237 Tag gradientTag;
238};
239
240enum DataTags { RHS = 0, LHS };
241
242struct MFrontInterface::CommonData::BlockData {
243
244 BlockData()
245 : isFiniteStrain(false), behaviourPath(DEFAULT_MFRONT_BEHAVIOUR_PATH),
246 behaviourName("LinearElasticity") {
247 dIssipation = 0;
248 storedEnergy = 0;
249 externalVariable = 0;
250 }
251
252 inline MoFEMErrorCode setTag(DataTags tag) {
253 MoFEMFunctionBeginHot;
254 if (tag == RHS) {
255 behDataPtr->K[0] = 0;
256 } else {
257 behDataPtr->K[0] = 5;
258 }
259 MoFEMFunctionReturnHot(0);
260 }
261
262 MoFEMErrorCode setBlockBehaviourData(bool set_params_from_blocks);
263
264 int iD;
265
266 bool isFiniteStrain;
267 string behaviourPath;
268 string behaviourName;
269
270 boost::shared_ptr<mgis::behaviour::Behaviour> mGisBehaviour;
271 mgis::behaviour::BehaviourDataView bView;
272 boost::shared_ptr<mgis::behaviour::BehaviourData> behDataPtr;
273
274 int sizeIntVar;
275 int sizeExtVar;
276 int sizeGradVar;
277 int sizeStressVar;
278
279 vector<double> params;
280
281 double dIssipation;
282 double storedEnergy;
283 double externalVariable;
284
285 Range eNts;
286};
287
288MoFEMErrorCode MFrontInterface::CommonData::BlockData::setBlockBehaviourData(
289 bool set_params_from_blocks) {
290 MoFEMFunctionBegin;
291 if (mGisBehaviour) {
292
293 auto &mgis_bv = *mGisBehaviour;
294
295 sizeIntVar = getArraySize(mgis_bv.isvs, mgis_bv.hypothesis);
296 sizeExtVar = getArraySize(mgis_bv.esvs, mgis_bv.hypothesis);
297 sizeGradVar = getArraySize(mgis_bv.gradients, mgis_bv.hypothesis);
298 sizeStressVar =
299 getArraySize(mgis_bv.thermodynamic_forces, mgis_bv.hypothesis);
300
301 behDataPtr = boost::make_shared<BehaviourData>(BehaviourData{mgis_bv});
302 bView = make_view(*behDataPtr);
303 const int total_number_of_params = mgis_bv.mps.size();
304
305 if (set_params_from_blocks) {
306
307 if (params.size() < total_number_of_params)
308 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
309 "Not enough parameters supplied for this block. We have %zu "
310 "provided where %d are necessary for this block",
311 params.size(), total_number_of_params);
312
313 for (int dd = 0; dd < total_number_of_params; ++dd) {
314 setMaterialProperty(behDataPtr->s0, dd, params[dd]);
315 setMaterialProperty(behDataPtr->s1, dd, params[dd]);
316 }
317 }
318
319 if (isFiniteStrain) {
320 behDataPtr->K[0] = 0; // no tangent
321 behDataPtr->K[1] = 2; // PK1
322 behDataPtr->K[2] = 2; // PK1
323 } else {
324 behDataPtr->K[0] = 0; // no tangent
325 behDataPtr->K[1] = 0; // cauchy
326 }
327
328 for (auto &mb : {&behDataPtr->s0, &behDataPtr->s1}) {
329 mb->dissipated_energy = dIssipation;
330 mb->stored_energy = storedEnergy;
331 setExternalStateVariable(*mb, 0, externalVariable);
332 }
333 }
334
335 MoFEMFunctionReturn(0);
336}
337
338MoFEMErrorCode MFrontInterface::CommonData::setBlocks(int dim) {
339 MoFEMFunctionBegin;
340 string block_name = "MAT_MFRONT";
341 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
342 if (it->getName().compare(0, block_name.size(), block_name) == 0) {
343 std::vector<double> block_data;
344 // FIXME: TODO: maybe this should be set only from the command line!!!
345 CHKERR it->getAttributes(block_data);
346 const int id = it->getMeshsetId();
347 EntityHandle meshset = it->getMeshset();
348 CHKERR mField.get_moab().get_entities_by_dimension(
349 meshset, dim, setOfBlocksData[id].eNts, true);
350 for (auto ent : setOfBlocksData[id].eNts)
351 blocksIDmap[ent] = id;
352
353 setOfBlocksData[id].iD = id;
354 setOfBlocksData[id].params.resize(block_data.size());
355
356 for (int n = 0; n != block_data.size(); n++)
357 setOfBlocksData[id].params[n] = block_data[n];
358 }
359 }
360
361 MoFEMFunctionReturn(0);
362}
363
364MoFEMErrorCode MFrontInterface::CommonData::getInternalVar(
365 const EntityHandle fe_ent, const int nb_gauss_pts, const int var_size,
366 const int grad_size, const int stress_size, bool is_large_strain) {
367 MoFEMFunctionBegin;
368
369 auto mget_tag_data = [&](Tag &m_tag, boost::shared_ptr<MatrixDouble> &m_mat,
370 const int &m_size, bool is_def_grad = false) {
371 MoFEMFunctionBeginHot;
372
373 double *tag_data;
374 int tag_size;
375 rval = mField.get_moab().tag_get_by_ptr(
376 m_tag, &fe_ent, 1, (const void **)&tag_data, &tag_size);
377
378 if (rval != MB_SUCCESS || tag_size != m_size * nb_gauss_pts) {
379 m_mat->resize(nb_gauss_pts, m_size, false);
380 m_mat->clear();
381 // initialize deformation gradient properly
382 if (is_def_grad && is_large_strain)
383 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
384 (*m_mat)(gg, 0) = 1;
385 (*m_mat)(gg, 1) = 1;
386 (*m_mat)(gg, 2) = 1;
387 }
388 void const *tag_data2[] = {&*m_mat->data().begin()};
389 const int tag_size2 = m_mat->data().size();
390 CHKERR mField.get_moab().tag_set_by_ptr(m_tag, &fe_ent, 1, tag_data2,
391 &tag_size2);
392 } else {
393 MatrixAdaptor tag_vec = MatrixAdaptor(
394 nb_gauss_pts, m_size,
395 ublas::shallow_array_adaptor<double>(tag_size, tag_data));
396
397 *m_mat = tag_vec;
398 }
399
400 MoFEMFunctionReturnHot(0);
401 };
402
403 CHKERR mget_tag_data(internalVariableTag, internalVariablePtr, var_size);
404 CHKERR mget_tag_data(stressTag, mPrevStressPtr, stress_size);
405 CHKERR mget_tag_data(gradientTag, mPrevGradPtr, grad_size, true);
406
407 MoFEMFunctionReturn(0);
408}
409
410MoFEMErrorCode
411MFrontInterface::CommonData::setInternalVar(const EntityHandle fe_ent) {
412 MoFEMFunctionBegin;
413
414 auto mset_tag_data = [&](Tag &m_tag, boost::shared_ptr<MatrixDouble> &m_mat) {
415 MoFEMFunctionBeginHot;
416 void const *tag_data[] = {&*m_mat->data().begin()};
417 const int tag_size = m_mat->data().size();
418 CHKERR mField.get_moab().tag_set_by_ptr(m_tag, &fe_ent, 1, tag_data,
419 &tag_size);
420 MoFEMFunctionReturnHot(0);
421 };
422
423 CHKERR mset_tag_data(internalVariableTag, internalVariablePtr);
424 CHKERR mset_tag_data(stressTag, mPrevStressPtr);
425 CHKERR mset_tag_data(gradientTag, mPrevGradPtr);
426
427 MoFEMFunctionReturn(0);
428}
429
430MoFEMErrorCode MFrontInterface::CommonData::createTags() {
431 MoFEMFunctionBegin;
432 const int default_length = 0;
433 CHKERR mField.get_moab().tag_get_handle(
434 "_INTERNAL_VAR", default_length, MB_TYPE_DOUBLE, internalVariableTag,
435 MB_TAG_CREAT | MB_TAG_VARLEN | MB_TAG_SPARSE, PETSC_NULLPTR);
436 CHKERR mField.get_moab().tag_get_handle(
437 "_STRESS_TAG", default_length, MB_TYPE_DOUBLE, stressTag,
438 MB_TAG_CREAT | MB_TAG_VARLEN | MB_TAG_SPARSE, PETSC_NULLPTR);
439 CHKERR mField.get_moab().tag_get_handle(
440 "_GRAD_TAG", default_length, MB_TYPE_DOUBLE, gradientTag,
441 MB_TAG_CREAT | MB_TAG_VARLEN | MB_TAG_SPARSE, PETSC_NULLPTR);
442
443 MoFEMFunctionReturn(0);
444}
445
446MoFEMErrorCode MFrontInterface::CommonData::clearTags() {
447 MoFEMFunctionBegin;
448 double zero = 0;
449 for (auto &[id, data] : setOfBlocksData) {
450 CHKERR mField.get_moab().tag_clear_data(internalVariableTag, data.eNts,
451 &zero);
452 CHKERR mField.get_moab().tag_clear_data(stressTag, data.eNts, &zero);
453 CHKERR mField.get_moab().tag_clear_data(gradientTag, data.eNts, &zero);
454 }
455 MoFEMFunctionReturn(0);
456}
457
458template <ModelHypothesis MH, AssemblyType AT>
459MFrontInterfaceImpl<MH, AT>::MFrontInterfaceImpl(MoFEM::Interface &m_field)
460 : mField(m_field) {
461
462 isFiniteKinematics = false;
463 saveGauss = PETSC_FALSE;
464 saveDomain = PETSC_TRUE;
465 testJacobian = PETSC_FALSE;
466 randomFieldScale = 1.0;
467 optionsPrefix = "mf_";
468 monitorPtr = nullptr;
469 fieldBase = NOBASE;
470
471 if (!LogManager::checkIfChannelExist("MFrontInterfaceWorld")) {
472 auto core_log = logging::core::get();
473
474 core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(),
475 "MFrontInterfaceWorld"));
476 core_log->add_sink(LogManager::createSink(LogManager::getStrmSync(),
477 "MFrontInterfaceSync"));
478 core_log->add_sink(LogManager::createSink(LogManager::getStrmSelf(),
479 "MFrontInterfaceSelf"));
480
481 LogManager::setLog("MFrontInterfaceWorld");
482 LogManager::setLog("MFrontInterfaceSync");
483 LogManager::setLog("MFrontInterfaceSelf");
484
485 MOFEM_LOG_TAG("MFrontInterfaceWorld", "MFrontInterface");
486 MOFEM_LOG_TAG("MFrontInterfaceSync", "MFrontInterface");
487 MOFEM_LOG_TAG("MFrontInterfaceSelf", "MFrontInterface");
488 }
489
490 MOFEM_LOG("MFrontInterfaceWorld", Sev::noisy) << "MFront Interface created";
491}
492
493template <ModelHypothesis MH, AssemblyType AT>
494MoFEMErrorCode MFrontInterfaceImpl<MH, AT>::getCommandLineParameters() {
495 MoFEMFunctionBegin;
496
497 PetscOptionsBegin(PETSC_COMM_WORLD, optionsPrefix.c_str(), "", "none");
498
499 CHKERR PetscOptionsBool("-save_gauss", "save gauss pts (internal variables)",
500 "", saveGauss, &saveGauss, PETSC_NULLPTR);
501 CHKERR PetscOptionsBool("-save_domain", "save results on a domain mesh", "",
502 saveDomain, &saveDomain, PETSC_NULLPTR);
503
504 CHKERR PetscOptionsBool("-test_jacobian", "test Jacobian (LHS matrix)", "",
505 testJacobian, &testJacobian, PETSC_NULLPTR);
506 CHKERR PetscOptionsReal("-random_field_scale",
507 "scale for the finite difference jacobian", "",
508 randomFieldScale, &randomFieldScale, PETSC_NULLPTR);
509
510 if (saveGauss)
511 moabGaussIntPtr = boost::shared_ptr<moab::Interface>(new moab::Core());
512
513 commonDataPtr = boost::make_shared<MFrontInterface::CommonData>(mField);
514
515 commonDataPtr->setBlocks(DIM);
516 commonDataPtr->createTags();
517
518 commonDataPtr->mGradPtr = boost::make_shared<MatrixDouble>();
519 commonDataPtr->mStressPtr = boost::make_shared<MatrixDouble>();
520 commonDataPtr->mFullStrainPtr = boost::make_shared<MatrixDouble>();
521 commonDataPtr->mFullStressPtr = boost::make_shared<MatrixDouble>();
522 commonDataPtr->mDispPtr = boost::make_shared<MatrixDouble>();
523 commonDataPtr->mPrevGradPtr = boost::make_shared<MatrixDouble>();
524 commonDataPtr->mPrevStressPtr = boost::make_shared<MatrixDouble>();
525 commonDataPtr->materialTangentPtr = boost::make_shared<MatrixDouble>();
526 commonDataPtr->mFullTangentPtr = boost::make_shared<MatrixDouble>();
527 commonDataPtr->internalVariablePtr = boost::make_shared<MatrixDouble>();
528
529 if (commonDataPtr->setOfBlocksData.empty() ||
530 commonDataPtr->blocksIDmap.empty()) {
531 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND,
532 "No blocksets on the mesh have been provided for MFront (e.g. "
533 "MAT_MFRONT_1)");
534 }
535
536 auto check_lib_finite_strain = [&](const std::string &lib,
537 const std::string &beh_name, bool &flag) {
538 MoFEMFunctionBeginHot;
539
540 ifstream f(lib.c_str());
541 if (!f.good())
542 MOFEM_LOG("MFrontInterfaceWorld", Sev::error)
543 << "Problem with the behaviour path: " << lib;
544
545 auto &lm = mgis::LibrariesManager::get();
546 flag = bool(lm.getBehaviourType(lib, beh_name) == 2) &&
547 (lm.getBehaviourKinematic(lib, beh_name) == 3);
548 MoFEMFunctionReturnHot(0);
549 };
550
551 auto op = FiniteStrainBehaviourOptions{};
552 op.stress_measure = FiniteStrainBehaviourOptions::PK1;
553 op.tangent_operator = FiniteStrainBehaviourOptions::DPK1_DF;
554
555 for (auto &block : commonDataPtr->setOfBlocksData) {
556 const int &id = block.first;
557 auto &lib_path = block.second.behaviourPath;
558 auto &name = block.second.behaviourName;
559 const string param_name = "-block_" + to_string(id);
560 const string param_path = "-lib_path_" + to_string(id);
561 const string param_from_blocks = "-params_" + to_string(id);
562 PetscBool set_from_blocks = PETSC_FALSE;
563 char char_name[255];
564 PetscBool is_param;
565
566 CHKERR PetscOptionsBool(param_from_blocks.c_str(),
567 "set parameters from blocks", "", set_from_blocks,
568 &set_from_blocks, PETSC_NULLPTR);
569
570 CHKERR PetscOptionsString(param_name.c_str(), "name of the behaviour", "",
571 "LinearElasticity", char_name, 255, &is_param);
572 if (is_param)
573 name = string(char_name);
574 string default_lib_path = DEFAULT_MFRONT_BEHAVIOUR_PATH;
575 CHKERR PetscOptionsString(
576 param_path.c_str(), "path to the behaviour library", "",
577 default_lib_path.c_str(), char_name, 255, &is_param);
578 if (is_param)
579 lib_path = string(char_name);
580 auto &mgis_bv_ptr = block.second.mGisBehaviour;
581 bool is_finite_strain = false;
582
583 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform)
584 << "Loading behaviour from " << lib_path;
585
586 CHKERR check_lib_finite_strain(lib_path, name, is_finite_strain);
587
588 mgis::behaviour::Hypothesis h;
589 switch (MH) {
590 case TRIDIMENSIONAL:
591 h = mgis::behaviour::Hypothesis::TRIDIMENSIONAL;
592 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform)
593 << "Model hypothesis: TRIDIMENSIONAL";
594 break;
595 case PLANESTRAIN:
596 h = mgis::behaviour::Hypothesis::PLANESTRAIN;
597 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform)
598 << "Model hypothesis: PLANESTRAIN";
599 break;
600 case AXISYMMETRICAL:
601 h = mgis::behaviour::Hypothesis::AXISYMMETRICAL;
602 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform)
603 << "Model hypothesis: AXISYMMETRICAL";
604 break;
605 default:
606 break;
607 }
608
609 if (is_finite_strain) {
610 mgis_bv_ptr = boost::make_shared<Behaviour>(load(op, lib_path, name, h));
611 block.second.isFiniteStrain = true;
612 } else
613 mgis_bv_ptr = boost::make_shared<Behaviour>(load(lib_path, name, h));
614
615 CHKERR block.second.setBlockBehaviourData(set_from_blocks);
616 for (size_t dd = 0; dd < mgis_bv_ptr->mps.size(); ++dd) {
617 double my_param = 0;
618 PetscBool is_set = PETSC_FALSE;
619 string param_cmd = "-param_" + to_string(id) + "_" + to_string(dd);
620 CHKERR PetscOptionsScalar(param_cmd.c_str(), "parameter from cmd", "",
621 my_param, &my_param, &is_set);
622 if (!is_set)
623 continue;
624 setMaterialProperty(block.second.behDataPtr->s0, dd, my_param);
625 setMaterialProperty(block.second.behDataPtr->s1, dd, my_param);
626 }
627
628 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform)
629 << mgis_bv_ptr->behaviour << " behaviour loaded on block "
630 << block.first;
631
632 if (is_finite_strain)
633 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform)
634 << "Finite Strain Kinematics";
635 else
636 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform)
637 << "Small Strain Kinematics";
638
639 if (mgis_bv_ptr->isvs.size()) {
640 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform) << "Internal variables:";
641 for (const auto &is : mgis_bv_ptr->isvs) {
642 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform) << ": " << is.name;
643 }
644 }
645
646 if (mgis_bv_ptr->esvs.size()) {
647 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform) << "External variables:";
648 for (const auto &es : mgis_bv_ptr->esvs) {
649 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform) << ": " << es.name;
650 }
651 }
652
653 auto it = block.second.behDataPtr->s0.material_properties.begin();
654 int nb = 0;
655
656 if (mgis_bv_ptr->mps.size()) {
657 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform) << "Material properties:";
658 for (const auto &mp : mgis_bv_ptr->mps) {
659 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform)
660 << nb++ << " : " << mp.name << " = " << *it++;
661 }
662 }
663
664 if (mgis_bv_ptr->params.size()) {
665 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform) << "Real parameters:";
666 for (const auto &p : mgis_bv_ptr->params) {
667 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform) << nb++ << " : " << p;
668 }
669 }
670
671 if (mgis_bv_ptr->iparams.size()) {
672 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform) << "Integer parameters:";
673 for (const auto &p : mgis_bv_ptr->iparams) {
674 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform) << nb++ << " : " << p;
675 }
676 }
677
678 if (mgis_bv_ptr->usparams.size()) {
679 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform)
680 << "Unsigned short parameters:";
681 for (const auto &p : mgis_bv_ptr->usparams) {
682 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform) << nb++ << " : " << p;
683 }
684 }
685 }
686
687 PetscOptionsEnd();
688
689 auto check_behaviours_kinematics = [&](bool &is_finite_kin) {
690 MoFEMFunctionBeginHot;
691 is_finite_kin =
692 commonDataPtr->setOfBlocksData.begin()->second.isFiniteStrain;
693 for (auto &block : commonDataPtr->setOfBlocksData) {
694 if (block.second.isFiniteStrain != is_finite_kin)
695 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
696 "All used MFront behaviours have to be of same kinematics "
697 "(small or "
698 "large strains)");
699 }
700 MoFEMFunctionReturnHot(0);
701 };
702
703 CHKERR check_behaviours_kinematics(isFiniteKinematics);
704
705 auto get_base = [&]() {
706 Range domain_ents;
707 CHKERR mField.get_moab().get_entities_by_dimension(0, DIM, domain_ents,
708 true);
709 if (domain_ents.empty())
710 CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Empty mesh");
711 const auto type = type_from_handle(domain_ents[0]);
712 switch (type) {
713 case MBQUAD:
714 return DEMKOWICZ_JACOBI_BASE;
715 case MBHEX:
716 return DEMKOWICZ_JACOBI_BASE;
717 case MBTRI:
718 return AINSWORTH_LEGENDRE_BASE;
719 case MBTET:
720 return AINSWORTH_LEGENDRE_BASE;
721 default:
722 CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Element type not handled");
723 }
724 return NOBASE;
725 };
726
727 fieldBase = get_base();
728 MOFEM_LOG("MFrontInterfaceWorld", Sev::verbose)
729 << "MFront post process projection base: "
730 << ApproximationBaseNames[fieldBase];
731
732 MoFEMFunctionReturn(0);
733}
734
735template <int DIM>
736using Tensor4Pack =
737 FTensor::Tensor4<FTensor::CursorPtr<double *, DIM * DIM * DIM * DIM>, DIM,
738 DIM, DIM, DIM>;
739
740template <int DIM>
741using Tensor2Pack =
742 FTensor::Tensor2<FTensor::CursorPtr<double *, DIM * DIM>, DIM, DIM>;
743
744template <int DIM>
745using Tensor1Pack = FTensor::Tensor1<FTensor::CursorPtr<double *, DIM>, DIM>;
746
747using Tensor1PackCoords = FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>;
748
749template <int DIM>
750using DdgPack = FTensor::Ddg<
751 FTensor::CursorPtr<double *,
752 ((DIM * (DIM + 1)) / 2) * ((DIM * (DIM + 1)) / 2)>,
753 DIM, DIM>;
754
755template <typename T> inline size_t get_paraview_size(T &vsize) {
756 return vsize > 1 ? (vsize > 3 ? 9 : 3) : 1;
757};
758
759
760template <typename T1, typename T2>
761inline MoFEMErrorCode get_tensor4_from_voigt(const T1 &K, T2 &D) {
762 MoFEMFunctionBeginHot;
763
764 Number<0> N0;
765 Number<1> N1;
766 Number<2> N2;
767
768 // 3D finite strain
769 if (std::is_same<T2, Tensor4Pack<3>>::value) {
770 D(N0, N0, N0, N0) = K[0];
771 D(N0, N0, N1, N1) = K[1];
772 D(N0, N0, N2, N2) = K[2];
773 D(N0, N0, N0, N1) = K[3];
774 D(N0, N0, N1, N0) = K[4];
775 D(N0, N0, N0, N2) = K[5];
776 D(N0, N0, N2, N0) = K[6];
777 D(N0, N0, N1, N2) = K[7];
778 D(N0, N0, N2, N1) = K[8];
779 D(N1, N1, N0, N0) = K[9];
780 D(N1, N1, N1, N1) = K[10];
781 D(N1, N1, N2, N2) = K[11];
782 D(N1, N1, N0, N1) = K[12];
783 D(N1, N1, N1, N0) = K[13];
784 D(N1, N1, N0, N2) = K[14];
785 D(N1, N1, N2, N0) = K[15];
786 D(N1, N1, N1, N2) = K[16];
787 D(N1, N1, N2, N1) = K[17];
788 D(N2, N2, N0, N0) = K[18];
789 D(N2, N2, N1, N1) = K[19];
790 D(N2, N2, N2, N2) = K[20];
791 D(N2, N2, N0, N1) = K[21];
792 D(N2, N2, N1, N0) = K[22];
793 D(N2, N2, N0, N2) = K[23];
794 D(N2, N2, N2, N0) = K[24];
795 D(N2, N2, N1, N2) = K[25];
796 D(N2, N2, N2, N1) = K[26];
797 D(N0, N1, N0, N0) = K[27];
798 D(N0, N1, N1, N1) = K[28];
799 D(N0, N1, N2, N2) = K[29];
800 D(N0, N1, N0, N1) = K[30];
801 D(N0, N1, N1, N0) = K[31];
802 D(N0, N1, N0, N2) = K[32];
803 D(N0, N1, N2, N0) = K[33];
804 D(N0, N1, N1, N2) = K[34];
805 D(N0, N1, N2, N1) = K[35];
806 D(N1, N0, N0, N0) = K[36];
807 D(N1, N0, N1, N1) = K[37];
808 D(N1, N0, N2, N2) = K[38];
809 D(N1, N0, N0, N1) = K[39];
810 D(N1, N0, N1, N0) = K[40];
811 D(N1, N0, N0, N2) = K[41];
812 D(N1, N0, N2, N0) = K[42];
813 D(N1, N0, N1, N2) = K[43];
814 D(N1, N0, N2, N1) = K[44];
815 D(N0, N2, N0, N0) = K[45];
816 D(N0, N2, N1, N1) = K[46];
817 D(N0, N2, N2, N2) = K[47];
818 D(N0, N2, N0, N1) = K[48];
819 D(N0, N2, N1, N0) = K[49];
820 D(N0, N2, N0, N2) = K[50];
821 D(N0, N2, N2, N0) = K[51];
822 D(N0, N2, N1, N2) = K[52];
823 D(N0, N2, N2, N1) = K[53];
824 D(N2, N0, N0, N0) = K[54];
825 D(N2, N0, N1, N1) = K[55];
826 D(N2, N0, N2, N2) = K[56];
827 D(N2, N0, N0, N1) = K[57];
828 D(N2, N0, N1, N0) = K[58];
829 D(N2, N0, N0, N2) = K[59];
830 D(N2, N0, N2, N0) = K[60];
831 D(N2, N0, N1, N2) = K[61];
832 D(N2, N0, N2, N1) = K[62];
833 D(N1, N2, N0, N0) = K[63];
834 D(N1, N2, N1, N1) = K[64];
835 D(N1, N2, N2, N2) = K[65];
836 D(N1, N2, N0, N1) = K[66];
837 D(N1, N2, N1, N0) = K[67];
838 D(N1, N2, N0, N2) = K[68];
839 D(N1, N2, N2, N0) = K[69];
840 D(N1, N2, N1, N2) = K[70];
841 D(N1, N2, N2, N1) = K[71];
842 D(N2, N1, N0, N0) = K[72];
843 D(N2, N1, N1, N1) = K[73];
844 D(N2, N1, N2, N2) = K[74];
845 D(N2, N1, N0, N1) = K[75];
846 D(N2, N1, N1, N0) = K[76];
847 D(N2, N1, N0, N2) = K[77];
848 D(N2, N1, N2, N0) = K[78];
849 D(N2, N1, N1, N2) = K[79];
850 D(N2, N1, N2, N1) = K[80];
851 }
852
853 // plane strain, finite strain
854 if (std::is_same<T2, Tensor4Pack<2>>::value) {
855 D(N0, N0, N0, N0) = K[0];
856 D(N0, N0, N1, N1) = K[1];
857 // K[2] is not used
858 D(N0, N0, N0, N1) = K[3];
859 D(N0, N0, N1, N0) = K[4];
860 D(N1, N1, N0, N0) = K[5];
861 D(N1, N1, N1, N1) = K[6];
862 // K[7] is not used
863 D(N1, N1, N0, N1) = K[8];
864 D(N1, N1, N1, N0) = K[9];
865 // K[10], K[11], K[12], K[13], K[14] are not used
866 D(N0, N1, N0, N0) = K[15];
867 D(N0, N1, N1, N1) = K[16];
868 // K[17] is not used
869 D(N0, N1, N0, N1) = K[18];
870 D(N0, N1, N1, N0) = K[19];
871 D(N1, N0, N0, N0) = K[20];
872 D(N1, N0, N1, N1) = K[21];
873 // K[22] is not used
874 D(N1, N0, N0, N1) = K[23];
875 D(N1, N0, N1, N0) = K[24];
876 }
877
878 // 3D small strain
879 if (std::is_same<T2, DdgPack<3>>::value) {
880 D(N0, N0, N0, N0) = K[0];
881 D(N0, N0, N1, N1) = K[1];
882 D(N0, N0, N2, N2) = K[2];
883
884 D(N0, N0, N0, N1) = inv_sqr2 * K[3];
885 D(N0, N0, N0, N2) = inv_sqr2 * K[4];
886 D(N0, N0, N1, N2) = inv_sqr2 * K[5];
887
888 D(N1, N1, N0, N0) = K[6];
889 D(N1, N1, N1, N1) = K[7];
890 D(N1, N1, N2, N2) = K[8];
891
892 D(N1, N1, N0, N1) = inv_sqr2 * K[9];
893 D(N1, N1, N0, N2) = inv_sqr2 * K[10];
894 D(N1, N1, N1, N2) = inv_sqr2 * K[11];
895
896 D(N2, N2, N0, N0) = K[12];
897 D(N2, N2, N1, N1) = K[13];
898 D(N2, N2, N2, N2) = K[14];
899
900 D(N2, N2, N0, N1) = inv_sqr2 * K[15];
901 D(N2, N2, N0, N2) = inv_sqr2 * K[16];
902 D(N2, N2, N1, N2) = inv_sqr2 * K[17];
903
904 D(N0, N1, N0, N0) = inv_sqr2 * K[18];
905 D(N0, N1, N1, N1) = inv_sqr2 * K[19];
906 D(N0, N1, N2, N2) = inv_sqr2 * K[20];
907
908 D(N0, N1, N0, N1) = 0.5 * K[21];
909 D(N0, N1, N0, N2) = 0.5 * K[22];
910 D(N0, N1, N1, N2) = 0.5 * K[23];
911
912 D(N0, N2, N0, N0) = inv_sqr2 * K[24];
913 D(N0, N2, N1, N1) = inv_sqr2 * K[25];
914 D(N0, N2, N2, N2) = inv_sqr2 * K[26];
915
916 D(N0, N2, N0, N1) = 0.5 * K[27];
917 D(N0, N2, N0, N2) = 0.5 * K[28];
918 D(N0, N2, N1, N2) = 0.5 * K[29];
919
920 D(N1, N2, N0, N0) = inv_sqr2 * K[30];
921 D(N1, N2, N1, N1) = inv_sqr2 * K[31];
922 D(N1, N2, N2, N2) = inv_sqr2 * K[32];
923
924 D(N1, N2, N0, N1) = 0.5 * K[33];
925 D(N1, N2, N0, N2) = 0.5 * K[34];
926 D(N1, N2, N1, N2) = 0.5 * K[35];
927 }
928
929 // plane strain, small strain
930 if (std::is_same<T2, DdgPack<2>>::value) {
931
932 D(N0, N0, N0, N0) = K[0];
933 D(N0, N0, N1, N1) = K[1];
934 // K[2] is not used
935
936 D(N0, N0, N0, N1) = inv_sqr2 * K[3];
937
938 D(N1, N1, N0, N0) = K[4];
939 D(N1, N1, N1, N1) = K[5];
940 // K[6] is not used
941
942 D(N1, N1, N0, N1) = inv_sqr2 * K[7];
943
944 // K[8], K[9], K[10], K[11] are not used
945
946 D(N0, N1, N0, N0) = inv_sqr2 * K[12];
947 D(N0, N1, N1, N1) = inv_sqr2 * K[13];
948
949 // K[14] is not used
950
951 D(N0, N1, N0, N1) = 0.5 * K[15];
952 }
953
954 MoFEMFunctionReturnHot(0);
955};
956
957template <bool IS_LARGE_STRAIN, typename T1, typename T2>
958inline MoFEMErrorCode get_full_tensor4_from_voigt(const T1 &K, T2 &D) {
959 MoFEMFunctionBeginHot;
960
961 Number<0> N0;
962 Number<1> N1;
963 Number<2> N2;
964
965 // 2D finite strain, full tensor
966 if constexpr (IS_LARGE_STRAIN) {
967 D(N0, N0, N0, N0) = K[0];
968 D(N0, N0, N1, N1) = K[1];
969 D(N0, N0, N2, N2) = K[2];
970 D(N0, N0, N0, N1) = K[3];
971 D(N0, N0, N1, N0) = K[4];
972 D(N1, N1, N0, N0) = K[5];
973 D(N1, N1, N1, N1) = K[6];
974 D(N1, N1, N2, N2) = K[7];
975 D(N1, N1, N0, N1) = K[8];
976 D(N1, N1, N1, N0) = K[9];
977 D(N2, N2, N0, N0) = K[10];
978 D(N2, N2, N1, N1) = K[11];
979 D(N2, N2, N2, N2) = K[12];
980 D(N2, N2, N0, N1) = K[13];
981 D(N2, N2, N1, N0) = K[14];
982 D(N0, N1, N0, N0) = K[15];
983 D(N0, N1, N1, N1) = K[16];
984 D(N0, N1, N2, N2) = K[17];
985 D(N0, N1, N0, N1) = K[18];
986 D(N0, N1, N1, N0) = K[19];
987 D(N1, N0, N0, N0) = K[20];
988 D(N1, N0, N1, N1) = K[21];
989 D(N1, N0, N2, N2) = K[22];
990 D(N1, N0, N0, N1) = K[23];
991 D(N1, N0, N1, N0) = K[24];
992 } else { // 2D small strain, full tensor
993
994 D(N0, N0, N0, N0) = K[0];
995 D(N0, N0, N1, N1) = K[1];
996 D(N0, N0, N2, N2) = K[2];
997
998 D(N0, N0, N0, N1) = inv_sqr2 * K[3];
999
1000 D(N1, N1, N0, N0) = K[4];
1001 D(N1, N1, N1, N1) = K[5];
1002 D(N1, N1, N2, N2) = K[6];
1003
1004 D(N1, N1, N0, N1) = inv_sqr2 * K[7];
1005
1006 D(N2, N2, N0, N0) = K[8];
1007 D(N2, N2, N1, N1) = K[9];
1008 D(N2, N2, N2, N2) = K[10];
1009
1010 D(N2, N2, N0, N1) = inv_sqr2 * K[11];
1011 D(N2, N2, N1, N0) = D(N2, N2, N0, N1);
1012
1013 D(N0, N1, N0, N0) = inv_sqr2 * K[12];
1014 D(N0, N1, N1, N1) = inv_sqr2 * K[13];
1015
1016 D(N0, N1, N2, N2) = inv_sqr2 * K[14];
1017 D(N1, N0, N2, N2) = D(N0, N1, N2, N2);
1018
1019 D(N0, N1, N0, N1) = 0.5 * K[15];
1020 }
1021
1022 MoFEMFunctionReturnHot(0);
1023};
1024
1025template <typename T> T get_tangent_tensor(MatrixDouble &mat);
1026
1027template <>
1028Tensor4Pack<3> get_tangent_tensor<Tensor4Pack<3>>(MatrixDouble &mat) {
1029 return getFTensor4FromMat<3, 3, 3, 3>(mat);
1030}
1031
1032template <>
1033Tensor4Pack<2> get_tangent_tensor<Tensor4Pack<2>>(MatrixDouble &mat) {
1034 return getFTensor4FromMat<2, 2, 2, 2>(mat);
1035}
1036
1037template <> DdgPack<3> get_tangent_tensor<DdgPack<3>>(MatrixDouble &mat) {
1038 return getFTensor4DdgFromMat<3, 3>(mat);
1039}
1040
1041template <> DdgPack<2> get_tangent_tensor<DdgPack<2>>(MatrixDouble &mat) {
1042 return getFTensor4DdgFromMat<2, 2>(mat);
1043}
1044
1045template <bool IS_LARGE_STRAIN, ModelHypothesis MH>
1046MoFEMErrorCode
1047mgis_integration(size_t gg, Tensor2Pack<MFrontEleType<MH>::SPACE_DIM> &t_grad,
1048 Tensor1Pack<MFrontEleType<MH>::SPACE_DIM> &t_disp,
1049 Tensor1PackCoords &t_coords,
1050 MFrontInterface::CommonData &common_data,
1051 MFrontInterface::CommonData::BlockData &block_data) {
1052 MoFEMFunctionBegin;
1053
1054 static constexpr int DIM = MFrontEleType<MH>::SPACE_DIM;
1055
1056 int check_integration;
1057 MatrixDouble &mat_int = *common_data.internalVariablePtr;
1058 MatrixDouble &mat_grad0 = *common_data.mPrevGradPtr;
1059 MatrixDouble &mat_stress0 = *common_data.mPrevStressPtr;
1060
1061 int &size_of_vars = block_data.sizeIntVar;
1062 int &size_of_grad = block_data.sizeGradVar;
1063 int &size_of_stress = block_data.sizeStressVar;
1064
1065 auto &mgis_bv = *block_data.mGisBehaviour;
1066
1067 FTensor::Index<'i', DIM> i;
1068 FTensor::Index<'j', DIM> j;
1069
1070 if constexpr (IS_LARGE_STRAIN) {
1071 FTensor::Tensor2<double, DIM, DIM> t_strain;
1072 t_strain(i, j) = t_grad(i, j) + kronecker_delta(i, j);
1073 if constexpr (MH == AXISYMMETRICAL) {
1074 setGradient(block_data.behDataPtr->s1, 0, size_of_grad,
1075 &*getVoigtVecAxisymm(t_strain, t_disp(0) / t_coords(0)).data());
1076 } else {
1077 setGradient(block_data.behDataPtr->s1, 0, size_of_grad,
1078 &*getVoigtVec<DIM>(t_strain).data());
1079 }
1080 } else {
1081 FTensor::Tensor2_symmetric<double, DIM> t_strain;
1082 t_strain(i, j) = (t_grad(i, j) || t_grad(j, i)) / 2;
1083 if constexpr (MH == AXISYMMETRICAL) {
1084 setGradient(
1085 block_data.behDataPtr->s1, 0, size_of_grad,
1086 &*getVoigtVecSymmAxisymm(t_strain, t_disp(0) / t_coords(0)).data());
1087 } else {
1088 setGradient(block_data.behDataPtr->s1, 0, size_of_grad,
1089 &*getVoigtVecSymm<DIM>(t_strain).data());
1090 }
1091 }
1092
1093 auto grad0_vec =
1094 getVectorAdaptor(&mat_grad0.data()[gg * size_of_grad], size_of_grad);
1095 setGradient(block_data.behDataPtr->s0, 0, size_of_grad, &*grad0_vec.begin());
1096
1097 auto stress0_vec = getVectorAdaptor(&mat_stress0.data()[gg * size_of_stress],
1098 size_of_stress);
1099 setThermodynamicForce(block_data.behDataPtr->s0, 0, size_of_stress,
1100 &*stress0_vec.begin());
1101
1102 if (size_of_vars) {
1103 auto internal_var =
1104 getVectorAdaptor(&mat_int.data()[gg * size_of_vars], size_of_vars);
1105 setInternalStateVariable(block_data.behDataPtr->s0, 0, size_of_vars,
1106 &*internal_var.begin());
1107 }
1108
1109 check_integration = mgis::behaviour::integrate(block_data.bView, mgis_bv);
1110 switch (check_integration) {
1111 case -1:
1112 MOFEM_LOG("WORLD", Sev::error) << "Mfront integration failed";
1113 break;
1114 case 0:
1115 MOFEM_LOG("WORLD", Sev::warning)
1116 << "Mfront integration succeeded but results are unreliable";
1117 break;
1118 case 1:
1119 default:
1120 break;
1121 }
1122
1123 MoFEMFunctionReturn(0);
1124}
1125
1126template <bool UPDATE, bool IS_LARGE_STRAIN, ModelHypothesis MH>
1127struct OpStressTmp : public MFrontEleType<MH>::DomainEleOp {
1128 static constexpr int DIM = MFrontEleType<MH>::SPACE_DIM;
1129 using DomainEleOp = typename MFrontEleType<MH>::DomainEleOp;
1130
1131 OpStressTmp(const std::string field_name,
1132 boost::shared_ptr<MFrontInterface::CommonData> common_data_ptr,
1133 boost::shared_ptr<FEMethod> monitor_ptr)
1134 : DomainEleOp(field_name, DomainEleOp::OPROW),
1135 commonDataPtr(common_data_ptr), monitorPtr(monitor_ptr) {
1136 std::fill(&DomainEleOp::doEntities[MBEDGE],
1137 &DomainEleOp::doEntities[MBMAXTYPE], false);
1138 }
1139 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
1140
1141private:
1142 boost::shared_ptr<MFrontInterface::CommonData> commonDataPtr;
1143 boost::shared_ptr<FEMethod> monitorPtr;
1144};
1145
1146template <ModelHypothesis MH>
1147using OpUpdateVariablesFiniteStrains = OpStressTmp<true, true, MH>;
1148
1149template <ModelHypothesis MH>
1150using OpUpdateVariablesSmallStrains = OpStressTmp<true, false, MH>;
1151
1152template <ModelHypothesis MH>
1153using OpStressFiniteStrains = OpStressTmp<false, true, MH>;
1154
1155template <ModelHypothesis MH>
1156using OpStressSmallStrains = OpStressTmp<false, false, MH>;
1157
1158template <bool UPDATE, bool IS_LARGE_STRAIN, ModelHypothesis MH>
1159MoFEMErrorCode OpStressTmp<UPDATE, IS_LARGE_STRAIN, MH>::doWork(int side,
1160 EntityType type,
1161 EntData &data) {
1162 MoFEMFunctionBegin;
1163
1164 Index<'i', 3> i;
1165 Index<'j', 3> j;
1166
1167 Index<'I', 2> I;
1168 Index<'J', 2> J;
1169
1170 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
1171 auto fe_ent = DomainEleOp::getNumeredEntFiniteElementPtr()->getEnt();
1172 auto id = commonDataPtr->blocksIDmap.at(fe_ent);
1173 auto &dAta = commonDataPtr->setOfBlocksData.at(id);
1174
1175 if (monitorPtr == nullptr)
1176 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_FOUND,
1177 "Time Monitor (FEMethod) has not been set for MFrontInterfaceImpl. "
1178 "Make sure to call setMonitor before calling "
1179 "opFactoryDomainRhs and opFactoryDomainLhs");
1180
1181 dAta.setTag(DataTags::RHS);
1182 dAta.behDataPtr->dt = monitorPtr->ts_dt;
1183 dAta.bView.dt = monitorPtr->ts_dt;
1184
1185 CHKERR commonDataPtr->getInternalVar(fe_ent, nb_gauss_pts, dAta.sizeIntVar,
1186 dAta.sizeGradVar, dAta.sizeStressVar,
1187 IS_LARGE_STRAIN);
1188
1189 MatrixDouble &mat_int = *commonDataPtr->internalVariablePtr;
1190 MatrixDouble &mat_grad0 = *commonDataPtr->mPrevGradPtr;
1191 MatrixDouble &mat_stress0 = *commonDataPtr->mPrevStressPtr;
1192
1193 auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->mGradPtr));
1194 auto t_disp = getFTensor1FromMat<DIM>(*(commonDataPtr->mDispPtr));
1195 auto t_coords = DomainEleOp::getFTensor1CoordsAtGaussPts();
1196
1197 commonDataPtr->mStressPtr->resize(nb_gauss_pts, DIM * DIM);
1198 commonDataPtr->mStressPtr->clear();
1199 auto t_stress = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->mStressPtr));
1200
1201 commonDataPtr->mFullStressPtr->resize(nb_gauss_pts, 3 * 3);
1202 commonDataPtr->mFullStressPtr->clear();
1203 auto t_full_stress =
1204 getFTensor2FromMat<3, 3>(*(commonDataPtr->mFullStressPtr));
1205
1206 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1207
1208 CHKERR mgis_integration<IS_LARGE_STRAIN, MH>(gg, t_grad, t_disp, t_coords,
1209 *commonDataPtr, dAta);
1210
1211 if constexpr (DIM == 3) {
1212 Tensor2<double, 3, 3> t_force;
1213 if constexpr (IS_LARGE_STRAIN) {
1214 t_force = Tensor2<double, 3, 3>(
1215 VOIGT_VEC_3D(getThermodynamicForce(dAta.behDataPtr->s1, 0)));
1216 } else {
1217 t_force = to_non_symm(Tensor2_symmetric<double, 3>(
1218 VOIGT_VEC_SYMM_3D(getThermodynamicForce(dAta.behDataPtr->s1, 0))));
1219 }
1220 t_stress(i, j) = t_force(i, j);
1221 } else if constexpr (DIM == 2) {
1222 Tensor2<double, 2, 2> t_force;
1223 Tensor2<double, 3, 3> t_full_force;
1224 if constexpr (IS_LARGE_STRAIN) {
1225 t_force = Tensor2<double, 2, 2>(
1226 VOIGT_VEC_2D(getThermodynamicForce(dAta.behDataPtr->s1, 0)));
1227 t_full_force = Tensor2<double, 3, 3>(
1228 VOIGT_VEC_2D_FULL(getThermodynamicForce(dAta.behDataPtr->s1, 0)));
1229 } else {
1230 t_force = to_non_symm(Tensor2_symmetric<double, 2>(
1231 VOIGT_VEC_SYMM_2D(getThermodynamicForce(dAta.behDataPtr->s1, 0))));
1232 t_full_force =
1233 to_non_symm(Tensor2_symmetric<double, 3>(VOIGT_VEC_SYMM_2D_FULL(
1234 getThermodynamicForce(dAta.behDataPtr->s1, 0))));
1235 }
1236 t_stress(I, J) = t_force(I, J);
1237 t_full_stress(i, j) = t_full_force(i, j);
1238 }
1239
1240 if constexpr (UPDATE) {
1241 for (int dd = 0; dd != dAta.sizeIntVar; ++dd) {
1242 mat_int(gg, dd) = *getInternalStateVariable(dAta.behDataPtr->s1, dd);
1243 }
1244 for (int dd = 0; dd != dAta.sizeGradVar; ++dd) {
1245 mat_grad0(gg, dd) = *getGradient(dAta.behDataPtr->s1, dd);
1246 }
1247 for (int dd = 0; dd != dAta.sizeStressVar; ++dd) {
1248 mat_stress0(gg, dd) = *getThermodynamicForce(dAta.behDataPtr->s1, dd);
1249 }
1250 }
1251
1252 ++t_stress;
1253 ++t_full_stress;
1254 ++t_grad;
1255 ++t_disp;
1256 ++t_coords;
1257 }
1258
1259 if constexpr (UPDATE) {
1260 CHKERR commonDataPtr->setInternalVar(fe_ent);
1261 }
1262
1263 MoFEMFunctionReturn(0);
1264}
1265
1266template <typename T, ModelHypothesis MH>
1267struct OpTangent : public MFrontEleType<MH>::DomainEleOp {
1268 static constexpr int DIM = MFrontEleType<MH>::SPACE_DIM;
1269 using DomainEleOp = typename MFrontEleType<MH>::DomainEleOp;
1270
1271 OpTangent(const std::string field_name,
1272 boost::shared_ptr<MFrontInterface::CommonData> common_data_ptr,
1273 boost::shared_ptr<FEMethod> monitor_ptr)
1274 : DomainEleOp(field_name, DomainEleOp::OPROW),
1275 commonDataPtr(common_data_ptr), monitorPtr(monitor_ptr) {
1276 std::fill(&DomainEleOp::doEntities[MBEDGE],
1277 &DomainEleOp::doEntities[MBMAXTYPE], false);
1278 }
1279 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
1280
1281private:
1282 boost::shared_ptr<MFrontInterface::CommonData> commonDataPtr;
1283 boost::shared_ptr<FEMethod> monitorPtr;
1284};
1285
1286template <ModelHypothesis MH>
1287using OpTangentFiniteStrains =
1288 struct OpTangent<Tensor4Pack<MFrontEleType<MH>::SPACE_DIM>, MH>;
1289
1290template <ModelHypothesis MH>
1291using OpTangentSmallStrains =
1292 struct OpTangent<DdgPack<MFrontEleType<MH>::SPACE_DIM>, MH>;
1293
1294template <typename T, ModelHypothesis MH>
1295MoFEMErrorCode OpTangent<T, MH>::doWork(int side, EntityType type,
1296 EntData &data) {
1297 MoFEMFunctionBegin;
1298
1299 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
1300 auto fe_ent = DomainEleOp::getNumeredEntFiniteElementPtr()->getEnt();
1301 auto id = commonDataPtr->blocksIDmap.at(fe_ent);
1302 auto &dAta = commonDataPtr->setOfBlocksData.at(id);
1303
1304 if (monitorPtr == nullptr)
1305 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_FOUND,
1306 "Time Monitor (FEMethod) has not been set for MFrontInterfaceImpl. "
1307 "Make sure to call setMonitor(monitor_ptr) before calling "
1308 "setOperators().");
1309
1310 dAta.setTag(DataTags::LHS);
1311 dAta.behDataPtr->dt = monitorPtr->ts_dt;
1312 dAta.bView.dt = monitorPtr->ts_dt;
1313
1314 constexpr bool IS_LARGE_STRAIN = std::is_same<T, Tensor4Pack<3>>::value ||
1315 std::is_same<T, Tensor4Pack<2>>::value;
1316
1317 CHKERR commonDataPtr->getInternalVar(fe_ent, nb_gauss_pts, dAta.sizeIntVar,
1318 dAta.sizeGradVar, dAta.sizeStressVar,
1319 IS_LARGE_STRAIN);
1320
1321 MatrixDouble &S_E = *(commonDataPtr->materialTangentPtr);
1322 MatrixDouble &F_E = *(commonDataPtr->mFullTangentPtr);
1323
1324 size_t tens_size = 36;
1325
1326 if constexpr (DIM == 2) {
1327 // plane strain
1328 if constexpr (IS_LARGE_STRAIN)
1329 tens_size = 16;
1330 else
1331 tens_size = 9;
1332 } else {
1333 if constexpr (IS_LARGE_STRAIN) // for finite strains
1334 tens_size = 81;
1335 }
1336
1337 S_E.resize(nb_gauss_pts, tens_size, false);
1338 auto D1 = get_tangent_tensor<T>(S_E);
1339
1340 size_t full_tens_size = 81;
1341 F_E.resize(nb_gauss_pts, full_tens_size, false);
1342 F_E.clear();
1343
1344 auto D2 = get_tangent_tensor<Tensor4Pack<3>>(F_E);
1345
1346 auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->mGradPtr));
1347 auto t_disp = getFTensor1FromMat<DIM>(*(commonDataPtr->mDispPtr));
1348 auto t_coords = DomainEleOp::getFTensor1CoordsAtGaussPts();
1349
1350 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1351
1352 CHKERR mgis_integration<IS_LARGE_STRAIN, MH>(gg, t_grad, t_disp, t_coords,
1353 *commonDataPtr, dAta);
1354
1355 CHKERR get_tensor4_from_voigt(&*dAta.behDataPtr->K.begin(), D1);
1356 CHKERR get_full_tensor4_from_voigt<IS_LARGE_STRAIN>(
1357 &*dAta.behDataPtr->K.begin(), D2);
1358
1359 ++D1;
1360 ++D2;
1361 ++t_grad;
1362 ++t_disp;
1363 ++t_coords;
1364 }
1365
1366 MoFEMFunctionReturn(0);
1367}
1368
1369template <AssemblyType AT>
1370struct OpAxisymmetricRhs
1371 : public FormsIntegrators<
1372 MFrontEleType<AXISYMMETRICAL>::DomainEleOp>::Assembly<AT>::OpBase {
1373 using DomainEleOp = typename MFrontEleType<AXISYMMETRICAL>::DomainEleOp;
1374 using OpBase = typename FormsIntegrators<DomainEleOp>::Assembly<AT>::OpBase;
1375
1376 OpAxisymmetricRhs(
1377 const std::string field_name,
1378 boost::shared_ptr<MFrontInterface::CommonData> common_data_ptr)
1379 : OpBase(field_name, field_name, DomainEleOp::OPROW),
1380 commonDataPtr(common_data_ptr) {};
1381
1382private:
1383 boost::shared_ptr<MFrontInterface::CommonData> commonDataPtr;
1384
1385 MoFEMErrorCode iNtegrate(EntData &row_data);
1386};
1387
1388template <AssemblyType AT>
1389MoFEMErrorCode OpAxisymmetricRhs<AT>::iNtegrate(EntData &row_data) {
1390 MoFEMFunctionBegin;
1391
1392 // get element volume
1393 const double vol = OpBase::getMeasure();
1394 // get integration weights
1395 auto t_w = OpBase::getFTensor0IntegrationWeight();
1396 // get coordinate at integration points
1397
1398 auto t_full_stress =
1399 getFTensor2FromMat<3, 3>(*(commonDataPtr->mFullStressPtr));
1400
1401 // loop over integration points
1402 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
1403
1404 auto t_nf = OpBase::template getNf<2>();
1405
1406 FTensor::Tensor0<double *> t_base(&row_data.getN()(gg, 0));
1407
1408 // take into account Jacobian
1409 const double alpha = t_w * vol * 2. * M_PI;
1410 // loop over rows base functions
1411 for (int rr = 0; rr != OpBase::nbRows / 2; ++rr) {
1412
1413 t_nf(0) += alpha * t_full_stress(2, 2) * t_base;
1414
1415 ++t_base;
1416 ++t_nf;
1417 }
1418
1419 ++t_full_stress;
1420 ++t_w; // move to another integration weight
1421 }
1422
1423 MoFEMFunctionReturn(0);
1424}
1425
1426template <AssemblyType AT>
1427struct OpAxisymmetricLhs
1428 : public FormsIntegrators<
1429 MFrontEleType<AXISYMMETRICAL>::DomainEleOp>::Assembly<AT>::OpBase {
1430 using DomainEleOp = typename MFrontEleType<AXISYMMETRICAL>::DomainEleOp;
1431 using OpBase = typename FormsIntegrators<DomainEleOp>::Assembly<AT>::OpBase;
1432
1433 OpAxisymmetricLhs(
1434 const std::string field_name,
1435 boost::shared_ptr<MFrontInterface::CommonData> common_data_ptr)
1436 : OpBase(field_name, field_name, DomainEleOp::OPROWCOL),
1437 commonDataPtr(common_data_ptr) {};
1438
1439private:
1440 boost::shared_ptr<MFrontInterface::CommonData> commonDataPtr;
1441
1442 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
1443};
1444
1445template <AssemblyType AT>
1446MoFEMErrorCode OpAxisymmetricLhs<AT>::iNtegrate(EntData &row_data,
1447 EntData &col_data) {
1448 MoFEMFunctionBegin;
1449
1450 // get element volume
1451 const double vol = OpBase::getMeasure();
1452 // get integration weights
1453 auto t_w = OpBase::getFTensor0IntegrationWeight();
1454 // get base function gradient on rows
1455 auto t_row_base = row_data.getFTensor0N();
1456 // get derivatives of base functions on rows
1457 auto t_row_diff_base = row_data.getFTensor1DiffN<2>();
1458 // get coordinate at integration points
1459 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
1460
1461 Number<0> N0;
1462 Number<1> N1;
1463 Number<2> N2;
1464
1465 auto t_D = getFTensor4FromMat<3, 3, 3, 3>(*(commonDataPtr->mFullTangentPtr));
1466
1467 // loop over integration points
1468 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
1469
1470 // Cylinder radius
1471 const double r_cylinder = t_coords(0);
1472
1473 // take into account Jacobean
1474 const double alpha = t_w * vol * 2. * M_PI;
1475
1476 // loop over rows base functions
1477 int rr = 0;
1478 for (; rr != OpBase::nbRows / 2; ++rr) {
1479
1480 // get sub matrix for the row
1481 auto t_m = OpBase::template getLocMat<2>(2 * rr);
1482
1483 // get derivatives of base functions for columns
1484 auto t_col_diff_base = col_data.getFTensor1DiffN<2>(gg, 0);
1485 // get base functions for columns
1486 auto t_col_base = col_data.getFTensor0N(gg, 0);
1487
1488 // loop over columns
1489 for (int cc = 0; cc != OpBase::nbCols / 2; cc++) {
1490
1491 t_m(0, 0) +=
1492 alpha * t_D(N0, N0, N2, N2) * t_col_base * t_row_diff_base(0);
1493
1494 t_m(1, 0) +=
1495 alpha * t_D(N1, N1, N2, N2) * t_col_base * t_row_diff_base(1);
1496
1497 t_m(0, 0) +=
1498 alpha * t_D(N2, N2, N0, N0) * t_col_diff_base(0) * t_row_base;
1499
1500 t_m(0, 1) +=
1501 alpha * t_D(N2, N2, N1, N1) * t_col_diff_base(1) * t_row_base;
1502
1503 t_m(0, 0) += alpha * (t_D(N2, N2, N2, N2) / r_cylinder) * t_col_base *
1504 t_row_base;
1505
1506 t_m(0, 0) +=
1507 alpha * t_D(N2, N2, N0, N1) * t_col_diff_base(1) * t_row_base;
1508
1509 t_m(0, 1) +=
1510 alpha * t_D(N2, N2, N1, N0) * t_col_diff_base(0) * t_row_base;
1511
1512 t_m(0, 0) +=
1513 alpha * t_D(N0, N1, N2, N2) * t_col_base * t_row_diff_base(1);
1514
1515 t_m(1, 0) +=
1516 alpha * t_D(N1, N0, N2, N2) * t_col_base * t_row_diff_base(0);
1517
1518 ++t_col_base;
1519 ++t_col_diff_base;
1520 ++t_m;
1521 }
1522
1523 ++t_row_base;
1524 ++t_row_diff_base;
1525 }
1526
1527 for (; rr < OpBase::nbRowBaseFunctions; ++rr) {
1528 ++t_row_base;
1529 ++t_row_diff_base;
1530 }
1531
1532 ++t_coords;
1533 ++t_w; // move to another integration weight
1534 ++t_D;
1535 }
1536
1537 MoFEMFunctionReturn(0);
1538}
1539
1540template <ModelHypothesis MH, AssemblyType AT>
1541MoFEMErrorCode MFrontInterfaceImpl<MH, AT>::opFactoryDomainRhs(
1542 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1543 std::string field_name) {
1544 MoFEMFunctionBegin;
1545
1546 auto jacobian = [&](const double r, const double, const double) {
1547 if (MH == AXISYMMETRICAL)
1548 return 2. * M_PI * r;
1549 else
1550 return 1.;
1551 };
1552
1553 auto add_domain_ops_rhs = [&](auto &pipeline) {
1554 if (isFiniteKinematics)
1555 pipeline.push_back(
1556 new OpStressFiniteStrains<MH>(field_name, commonDataPtr, monitorPtr));
1557 else
1558 pipeline.push_back(
1559 new OpStressSmallStrains<MH>(field_name, commonDataPtr, monitorPtr));
1560
1561 pipeline.push_back(
1562 new OpInternalForce(field_name, commonDataPtr->mStressPtr, jacobian));
1563
1564 if (MH == AXISYMMETRICAL)
1565 pipeline.push_back(new OpAxisymmetricRhs<AT>(field_name, commonDataPtr));
1566 };
1567
1568 auto add_domain_base_ops = [&](auto &pipeline) {
1569 pipeline.push_back(new OpCalculateVectorFieldValues<DIM>(
1570 field_name, commonDataPtr->mDispPtr));
1571 pipeline.push_back(new OpCalculateVectorFieldGradient<DIM, DIM>(
1572 field_name, commonDataPtr->mGradPtr));
1573 };
1574
1575 add_domain_base_ops(pip);
1576 add_domain_ops_rhs(pip);
1577
1578 MoFEMFunctionReturn(0);
1579}
1580
1581template <ModelHypothesis MH, AssemblyType AT>
1582MoFEMErrorCode MFrontInterfaceImpl<MH, AT>::opFactoryDomainLhs(
1583 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1584 std::string field_name) {
1585 MoFEMFunctionBegin;
1586
1587 auto jacobian = [&](const double r, const double, const double) {
1588 if (MH == AXISYMMETRICAL)
1589 return 2. * M_PI * r;
1590 else
1591 return 1.;
1592 };
1593
1594 auto add_domain_ops_lhs = [&](auto &pipeline) {
1595 if (isFiniteKinematics) {
1596 pipeline.push_back(new OpTangentFiniteStrains<MH>(
1597 field_name, commonDataPtr, monitorPtr));
1598 pipeline.push_back(new OpAssembleLhsFiniteStrains(
1599 field_name, field_name, commonDataPtr->materialTangentPtr, jacobian));
1600 } else {
1601 pipeline.push_back(
1602 new OpTangentSmallStrains<MH>(field_name, commonDataPtr, monitorPtr));
1603 pipeline.push_back(new OpAssembleLhsSmallStrains(
1604 field_name, field_name, commonDataPtr->materialTangentPtr, nullptr,
1605 jacobian));
1606 }
1607 if (MH == AXISYMMETRICAL)
1608 pipeline.push_back(new OpAxisymmetricLhs<AT>(field_name, commonDataPtr));
1609 };
1610
1611 auto add_domain_base_ops = [&](auto &pipeline) {
1612 pipeline.push_back(new OpCalculateVectorFieldValues<DIM>(
1613 field_name, commonDataPtr->mDispPtr));
1614 pipeline.push_back(new OpCalculateVectorFieldGradient<DIM, DIM>(
1615 field_name, commonDataPtr->mGradPtr));
1616 };
1617
1618 add_domain_base_ops(pip);
1619 add_domain_ops_lhs(pip);
1620
1621 MoFEMFunctionReturn(0);
1622}
1623
1624template <ModelHypothesis MH>
1625struct OpSaveGaussPts : public MFrontEleType<MH>::DomainEleOp {
1626 static constexpr int DIM = MFrontEleType<MH>::SPACE_DIM;
1627 using DomainEleOp = typename MFrontEleType<MH>::DomainEleOp;
1628
1629 OpSaveGaussPts(const std::string field_name, moab::Interface &moab_mesh,
1630 boost::shared_ptr<MFrontInterface::CommonData> common_data_ptr)
1631 : DomainEleOp(field_name, DomainEleOp::OPROW), internalVarMesh(moab_mesh),
1632 commonDataPtr(common_data_ptr), fieldName(field_name) {
1633 std::fill(&DomainEleOp::doEntities[MBEDGE],
1634 &DomainEleOp::doEntities[MBMAXTYPE], false);
1635 }
1636
1637 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1638 MoFEMFunctionBegin;
1639
1640 auto fe_ent = DomainEleOp::getNumeredEntFiniteElementPtr()->getEnt();
1641 auto id = commonDataPtr->blocksIDmap.at(fe_ent);
1642 auto &dAta = commonDataPtr->setOfBlocksData.at(id);
1643 auto &mgis_bv = *dAta.mGisBehaviour;
1644
1645 int &size_of_vars = dAta.sizeIntVar;
1646 int &size_of_grad = dAta.sizeGradVar;
1647 int &size_of_stress = dAta.sizeStressVar;
1648
1649 auto get_tag = [&](std::string name, size_t size) {
1650 std::array<double, 9> def;
1651 std::fill(def.begin(), def.end(), 0);
1652 Tag th;
1653 CHKERR internalVarMesh.tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE,
1654 th, MB_TAG_CREAT | MB_TAG_SPARSE,
1655 def.data());
1656 return th;
1657 };
1658
1659 auto t_stress = getFTensor2FromMat<3, 3>(*(commonDataPtr->mStressPtr));
1660
1661 MatrixDouble3by3 mat(3, 3);
1662
1663 auto th_disp = get_tag(fieldName, 3);
1664 auto th_stress = get_tag(mgis_bv.thermodynamic_forces[0].name, 9);
1665 auto th_grad = get_tag(mgis_bv.gradients[0].name, 9);
1666
1667 size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
1668 auto t_grad = getFTensor2FromMat<3, 3>(*(commonDataPtr->mGradPtr));
1669 auto t_disp = getFTensor1FromMat<3>(*(commonDataPtr->mDispPtr));
1670 CHKERR commonDataPtr->getInternalVar(fe_ent, nb_gauss_pts, size_of_vars,
1671 size_of_grad, size_of_stress);
1672
1673 MatrixDouble &mat_int = *commonDataPtr->internalVariablePtr;
1674 vector<Tag> tags_vec;
1675 bool is_large_strain = dAta.isFiniteStrain;
1676
1677 for (auto c : mgis_bv.isvs) {
1678 auto vsize = getVariableSize(c, mgis_bv.hypothesis);
1679 const size_t parav_siz = get_paraview_size(vsize);
1680 tags_vec.emplace_back(get_tag(c.name, parav_siz));
1681 }
1682
1683 if (!(side == 0 && type == MBVERTEX))
1684 MoFEMFunctionReturnHot(0);
1685
1686 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1687
1688 double coords[] = {0, 0, 0};
1689 EntityHandle vertex;
1690 for (int dd = 0; dd != 3; dd++)
1691 coords[dd] = DomainEleOp::getCoordsAtGaussPts()(gg, dd);
1692
1693 CHKERR internalVarMesh.create_vertex(coords, vertex);
1694 VectorDouble disps({t_disp(0), t_disp(1), t_disp(2)});
1695
1696 auto it = tags_vec.begin();
1697 for (auto c : mgis_bv.isvs) {
1698 auto vsize = getVariableSize(c, mgis_bv.hypothesis);
1699 const size_t parav_siz = get_paraview_size(vsize);
1700 const auto offset =
1701 getVariableOffset(mgis_bv.isvs, c.name, mgis_bv.hypothesis);
1702 auto vec =
1703 getVectorAdaptor(&mat_int.data()[gg * size_of_vars], size_of_vars);
1704 VectorDouble tag_vec = getVectorAdaptor(&vec[offset], vsize);
1705 tag_vec.resize(parav_siz);
1706
1707 CHKERR internalVarMesh.tag_set_data(*it, &vertex, 1, &*tag_vec.begin());
1708
1709 it++;
1710 }
1711
1712 // keep the convention consistent for postprocessing!
1713 array<double, 9> my_stress_vec{
1714 t_stress(0, 0), t_stress(1, 1), t_stress(2, 2),
1715 t_stress(0, 1), t_stress(1, 0), t_stress(0, 2),
1716 t_stress(2, 0), t_stress(1, 2), t_stress(2, 1)};
1717
1718 FTensor::Index<'i', 3> i;
1719 FTensor::Index<'j', 3> j;
1720
1721 array<double, 9> grad1_vec;
1722 if (is_large_strain) {
1723 FTensor::Tensor2<double, 3, 3> t_strain;
1724 t_strain(i, j) = t_grad(i, j) + kronecker_delta(i, j);
1725 grad1_vec = getVoigtVec<3>(t_strain);
1726 } else {
1727 FTensor::Tensor2_symmetric<double, 3> t_strain;
1728 t_strain(i, j) = (t_grad(i, j) || t_grad(j, i)) / 2;
1729 grad1_vec = getVoigtVecSymm<3>(t_strain);
1730 }
1731
1732 CHKERR internalVarMesh.tag_set_data(th_stress, &vertex, 1,
1733 my_stress_vec.data());
1734 CHKERR internalVarMesh.tag_set_data(th_grad, &vertex, 1,
1735 grad1_vec.data());
1736 CHKERR internalVarMesh.tag_set_data(th_disp, &vertex, 1, &*disps.begin());
1737
1738 ++t_grad;
1739 ++t_stress;
1740 ++t_disp;
1741 }
1742
1743 MoFEMFunctionReturn(0);
1744 }
1745
1746private:
1747 boost::shared_ptr<MFrontInterface::CommonData> commonDataPtr;
1748 moab::Interface &internalVarMesh;
1749 std::string fieldName;
1750};
1751
1752template <ModelHypothesis MH, AssemblyType AT>
1753MoFEMErrorCode MFrontInterfaceImpl<MH, AT>::setUpdateInternalVariablesOperators(
1754 ForcesAndSourcesCore::RuleHookFun rule, std::string field_name) {
1755 MoFEMFunctionBegin;
1756
1757 updateIntVariablesElePtr = boost::make_shared<DomainEle>(mField);
1758
1759 updateIntVariablesElePtr->getRuleHook = rule;
1760
1761 CHKERR AddHOOps<DIM, DIM, DIM>::add(
1762 updateIntVariablesElePtr->getOpPtrVector(), {H1});
1763
1764 updateIntVariablesElePtr->getOpPtrVector().push_back(
1765 new OpCalculateVectorFieldGradient<DIM, DIM>(field_name,
1766 commonDataPtr->mGradPtr));
1767 updateIntVariablesElePtr->getOpPtrVector().push_back(
1768 new OpCalculateVectorFieldValues<DIM>(field_name,
1769 commonDataPtr->mDispPtr));
1770 if (isFiniteKinematics)
1771 updateIntVariablesElePtr->getOpPtrVector().push_back(
1772 new OpUpdateVariablesFiniteStrains<MH>(field_name, commonDataPtr,
1773 monitorPtr));
1774 else
1775 updateIntVariablesElePtr->getOpPtrVector().push_back(
1776 new OpUpdateVariablesSmallStrains<MH>(field_name, commonDataPtr,
1777 monitorPtr));
1778 if (saveGauss && (MH == TRIDIMENSIONAL)) {
1779 auto &moab_gauss = *moabGaussIntPtr;
1780 updateIntVariablesElePtr->getOpPtrVector().push_back(
1781 new OpSaveGaussPts<MH>(field_name, moab_gauss, commonDataPtr));
1782 }
1783
1784 MoFEMFunctionReturn(0);
1785}
1786
1787template <bool IS_LARGE_STRAIN, ModelHypothesis MH>
1788struct OpSaveStress : public MFrontEleType<MH>::DomainEleOp {
1789 static constexpr int DIM = MFrontEleType<MH>::SPACE_DIM;
1790 using DomainEleOp = typename MFrontEleType<MH>::DomainEleOp;
1791
1792 OpSaveStress(const std::string field_name,
1793 boost::shared_ptr<MFrontInterface::CommonData> common_data_ptr)
1794 : DomainEleOp(field_name, DomainEleOp::OPROW),
1795 commonDataPtr(common_data_ptr) {
1796 std::fill(&DomainEleOp::doEntities[MBEDGE],
1797 &DomainEleOp::doEntities[MBMAXTYPE], false);
1798 }
1799 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1800 MoFEMFunctionBegin;
1801
1802 const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
1803 auto fe_ent = DomainEleOp::getNumeredEntFiniteElementPtr()->getEnt();
1804 auto id = commonDataPtr->blocksIDmap.at(fe_ent);
1805 auto &dAta = commonDataPtr->setOfBlocksData.at(id);
1806
1807 int &size_of_stress = dAta.sizeStressVar;
1808 int &size_of_grad = dAta.sizeGradVar;
1809
1810 CHKERR commonDataPtr->getInternalVar(fe_ent, nb_gauss_pts, dAta.sizeIntVar,
1811 dAta.sizeGradVar, dAta.sizeStressVar,
1812 IS_LARGE_STRAIN);
1813
1814 MatrixDouble &mat_stress = *commonDataPtr->mPrevStressPtr;
1815 MatrixDouble &mat_grad = *commonDataPtr->mPrevGradPtr;
1816
1817 commonDataPtr->mFullStressPtr->resize(nb_gauss_pts, 3 * 3);
1818 commonDataPtr->mFullStressPtr->clear();
1819 auto t_full_stress =
1820 getFTensor2FromMat<3, 3>(*(commonDataPtr->mFullStressPtr));
1821
1822 commonDataPtr->mFullStrainPtr->resize(nb_gauss_pts, 3 * 3);
1823 commonDataPtr->mFullStrainPtr->clear();
1824 auto t_full_strain =
1825 getFTensor2FromMat<3, 3>(*(commonDataPtr->mFullStrainPtr));
1826
1827 Index<'i', 3> i;
1828 Index<'j', 3> j;
1829
1830 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1831
1832 auto stress_vec = getVectorAdaptor(
1833 &mat_stress.data()[gg * size_of_stress], size_of_stress);
1834 auto grad_vec =
1835 getVectorAdaptor(&mat_grad.data()[gg * size_of_grad], size_of_grad);
1836
1837 Tensor2<double, 3, 3> t_stress;
1838 Tensor2<double, 3, 3> t_grad;
1839
1840 if constexpr (IS_LARGE_STRAIN) {
1841 if constexpr (DIM == 3) {
1842 t_stress = Tensor2<double, 3, 3>(VOIGT_VEC_3D(stress_vec));
1843 t_grad = Tensor2<double, 3, 3>(VOIGT_VEC_3D(grad_vec));
1844 } else if constexpr (DIM == 2) {
1845 t_stress = Tensor2<double, 3, 3>(VOIGT_VEC_2D_FULL(stress_vec));
1846 t_grad = Tensor2<double, 3, 3>(VOIGT_VEC_2D_FULL(grad_vec));
1847 }
1848 } else {
1849 if constexpr (DIM == 3) {
1850 t_stress = to_non_symm(
1851 Tensor2_symmetric<double, 3>(VOIGT_VEC_SYMM_3D(stress_vec)));
1852 t_grad = to_non_symm(
1853 Tensor2_symmetric<double, 3>(VOIGT_VEC_SYMM_3D(grad_vec)));
1854 } else if constexpr (DIM == 2) {
1855 t_stress = to_non_symm(
1856 Tensor2_symmetric<double, 3>(VOIGT_VEC_SYMM_2D_FULL(stress_vec)));
1857 t_grad = to_non_symm(
1858 Tensor2_symmetric<double, 3>(VOIGT_VEC_SYMM_2D_FULL(grad_vec)));
1859 }
1860 }
1861
1862 t_full_stress(i, j) = t_stress(i, j);
1863 t_full_strain(i, j) = t_grad(i, j);
1864
1865 ++t_full_stress;
1866 ++t_full_strain;
1867 }
1868
1869 MoFEMFunctionReturn(0);
1870 }
1871
1872private:
1873 boost::shared_ptr<MFrontInterface::CommonData> commonDataPtr;
1874};
1875
1876template <ModelHypothesis MH, AssemblyType AT>
1877MoFEMErrorCode MFrontInterfaceImpl<MH, AT>::setPostProcessOperators(
1878 ForcesAndSourcesCore::RuleHookFun rule, std::string fe_name,
1879 std::string field_name, int order) {
1880 MoFEMFunctionBegin;
1881
1882 postProcFe = boost::make_shared<PostProcDomainEle>(mField);
1883
1884 auto &pip = postProcFe->getOpPtrVector();
1885
1886 pip.push_back(new OpCalculateVectorFieldValues<DIM>(field_name,
1887 commonDataPtr->mDispPtr));
1888
1889 auto entity_data_l2 = boost::make_shared<EntitiesFieldData>(MBENTITYSET);
1890 auto mass_ptr = boost::make_shared<MatrixDouble>();
1891 auto strain_coeffs_ptr = boost::make_shared<MatrixDouble>();
1892 auto stress_coeffs_ptr = boost::make_shared<MatrixDouble>();
1893
1894 auto op_this = new OpLoopThis<DomainEle>(mField, fe_name, Sev::noisy);
1895 pip.push_back(op_this);
1896 pip.push_back(new OpDGProjectionEvaluation(
1897 commonDataPtr->mFullStrainPtr, strain_coeffs_ptr, entity_data_l2,
1898 fieldBase, L2, Sev::noisy,
1899 OpDGProjectionEvaluation::OutputLayout::GaussByCoeffs));
1900 pip.push_back(new OpDGProjectionEvaluation(
1901 commonDataPtr->mFullStressPtr, stress_coeffs_ptr, entity_data_l2,
1902 fieldBase, L2, Sev::noisy,
1903 OpDGProjectionEvaluation::OutputLayout::GaussByCoeffs));
1904
1905 auto fe_physics_ptr = op_this->getThisFEPtr();
1906 fe_physics_ptr->getRuleHook = rule;
1907
1908 fe_physics_ptr->getOpPtrVector().push_back(new OpDGProjectionMassMatrix(
1909 order, mass_ptr, entity_data_l2, fieldBase, L2));
1910 if (isFiniteKinematics) {
1911 fe_physics_ptr->getOpPtrVector().push_back(
1912 new OpSaveStress<true, MH>(field_name, commonDataPtr));
1913 } else {
1914 fe_physics_ptr->getOpPtrVector().push_back(
1915 new OpSaveStress<false, MH>(field_name, commonDataPtr));
1916 }
1917 fe_physics_ptr->getOpPtrVector().push_back(new OpDGProjectionCoefficients(
1918 commonDataPtr->mFullStrainPtr, strain_coeffs_ptr, mass_ptr,
1919 entity_data_l2, fieldBase, L2, Sev::noisy, DataLayout::GaussByCoeffs));
1920 fe_physics_ptr->getOpPtrVector().push_back(new OpDGProjectionCoefficients(
1921 commonDataPtr->mFullStressPtr, stress_coeffs_ptr, mass_ptr,
1922 entity_data_l2, fieldBase, L2, Sev::noisy, DataLayout::GaussByCoeffs));
1923
1924 using OpPPMapVec = OpPostProcMapInMoab<DIM, DIM>;
1925 using OpPPMapTen = OpPostProcMapInMoab<3, 3>;
1926
1927 pip.push_back(new OpPPMapVec(
1928
1929 postProcFe->getPostProcMesh(), postProcFe->getMapGaussPts(),
1930
1931 {},
1932
1933 {{field_name, commonDataPtr->mDispPtr}},
1934
1935 {},
1936
1937 {}
1938
1939 ));
1940
1941 pip.push_back(new OpPPMapTen(
1942
1943 postProcFe->getPostProcMesh(), postProcFe->getMapGaussPts(),
1944
1945 {},
1946
1947 {},
1948
1949 {{"STRAIN", commonDataPtr->mFullStrainPtr},
1950 {"STRESS", commonDataPtr->mFullStressPtr}},
1951
1952 {}
1953
1954 ));
1955
1956 MoFEMFunctionReturn(0);
1957}
1958
1959template <ModelHypothesis MH, AssemblyType AT>
1960MoFEMErrorCode MFrontInterfaceImpl<MH, AT>::postProcess(int step,
1961 SmartPetscObj<DM> dm,
1962 string fe_name) {
1963 MoFEMFunctionBegin;
1964
1965 auto make_vtks = [&]() {
1966 MoFEMFunctionBegin;
1967
1968 if (saveDomain) {
1969 CHKERR DMoFEMLoopFiniteElements(dm, fe_name, postProcFe);
1970 CHKERR postProcFe->writeFile("out_" + optionsPrefix +
1971 boost::lexical_cast<std::string>(step) +
1972 ".h5m");
1973 }
1974
1975 if (saveGauss) {
1976 string file_name = "out_" + optionsPrefix + "gauss_" +
1977 boost::lexical_cast<std::string>(step) + ".h5m";
1978
1979 CHKERR moabGaussIntPtr->write_file(file_name.c_str(), "MOAB",
1980 "PARALLEL=WRITE_PART");
1981 CHKERR moabGaussIntPtr->delete_mesh();
1982 }
1983
1984 MoFEMFunctionReturn(0);
1985 };
1986
1987 CHKERR make_vtks();
1988
1989 MoFEMFunctionReturn(0);
1990};
1991
1992template <ModelHypothesis MH, AssemblyType AT>
1993MoFEMErrorCode
1994MFrontInterfaceImpl<MH, AT>::updateInternalVariables(SmartPetscObj<DM> dm,
1995 std::string fe_name) {
1996 MoFEMFunctionBegin;
1997 CHKERR DMoFEMLoopFiniteElements(dm, fe_name, updateIntVariablesElePtr);
1998 MoFEMFunctionReturn(0);
1999}
2000
2001template struct MFrontInterfaceImpl<TRIDIMENSIONAL, AssemblyType::PETSC>;
2002template struct MFrontInterfaceImpl<AXISYMMETRICAL, AssemblyType::PETSC>;
2003template struct MFrontInterfaceImpl<PLANESTRAIN, AssemblyType::PETSC>;
2004
2005template struct MFrontInterfaceImpl<TRIDIMENSIONAL, AssemblyType::BLOCK_SCHUR>;
2006template struct MFrontInterfaceImpl<AXISYMMETRICAL, AssemblyType::BLOCK_SCHUR>;
2007template struct MFrontInterfaceImpl<PLANESTRAIN, AssemblyType::BLOCK_SCHUR>;
2008
2009} // namespace MoFEM
2010
2011#else // WITH_MGIS
2012
2013namespace MoFEM {
2014
2018
2019boost::shared_ptr<MFrontInterface>
2020createMFrontInterface(MoFEM::Interface &m_field, ModelHypothesis mh,
2021 AssemblyType at) {
2022 CHK_THROW_MESSAGE(
2023 MOFEM_NOT_IMPLEMENTED,
2024 "MoFEM built without MFront support. MFrontInterface is not available.");
2025}
2026
2027} // namespace MoFEM
2028
2029#endif // WITH_MGIS
setBlocks
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData > > &block_sets_ptr)
type
std::string type
Definition JsonConfigManager.cpp:28
MFrontInterface.hpp
MFrontInterface.
SPACE_DIM
constexpr int SPACE_DIM
Definition child_and_parent.cpp:17
DomainEleOp
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
Definition child_and_parent.cpp:37
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition child_and_parent.cpp:35
FTensor::Tensor2_symmetric
Definition Tensor2_symmetric_value.hpp:14
FieldApproximationBase
FieldApproximationBase
approximation base
Definition definitions.h:58
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
NOBASE
@ NOBASE
Definition definitions.h:59
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition definitions.h:612
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
L2
@ L2
field with C-1 continuity
Definition definitions.h:88
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
BLOCKSET
@ BLOCKSET
Definition definitions.h:161
MOFEM_NOT_FOUND
@ MOFEM_NOT_FOUND
Definition definitions.h:33
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
ApproximationBaseNames
static const char *const ApproximationBaseNames[]
Definition definitions.h:72
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:456
order
constexpr int order
Definition dg_projection.cpp:38
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:576
MoFEM::AssemblyType
AssemblyType
[Storage and set boundary conditions]
Definition FormsIntegrators.hpp:200
MoFEM::PETSC
@ PETSC
Standard PETSc assembly.
Definition FormsIntegrators.hpp:201
MoFEM::BLOCK_MAT
@ BLOCK_MAT
Block matrix assembly.
Definition FormsIntegrators.hpp:203
MoFEM::SCHUR
@ SCHUR
Schur complement assembly.
Definition FormsIntegrators.hpp:202
MoFEM::BLOCK_SCHUR
@ BLOCK_SCHUR
Block Schur assembly.
Definition FormsIntegrators.hpp:204
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:316
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition LogManager.hpp:347
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
Definition MeshsetsManager.hpp:94
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
c
const double c
speed of light (cm/ns)
Definition initial_diffusion.cpp:40
D
double D
Definition initial_diffusion.cpp:45
n
const double n
refractive index of diffusive medium
Definition initial_diffusion.cpp:39
J
FTensor::Index< 'J', DIM1 > J
Definition level_set.cpp:30
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:21
FTensor
Tensors class implemented by Walter Landry.
Definition FTensor.hpp:51
FTensor::kronecker_delta
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
Definition Kronecker_Delta.hpp:81
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition ddTensor0.hpp:33
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition Exceptions.hpp:74
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition Exceptions.hpp:56
MoFEM::Types::MatrixAdaptor
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Matrix adaptor.
Definition Types.hpp:131
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition Types.hpp:77
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition Types.hpp:68
MoFEM::Types::MatrixDouble3by3
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition Types.hpp:105
MoFEM
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
MoFEM::getVectorAdaptor
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition Templates.hpp:49
MoFEM::type_from_handle
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition Templates.hpp:2548
MoFEM::getVoigtVecAxisymm
auto getVoigtVecAxisymm(T &t_mat, const double hoop_term)
Definition Templates.hpp:1911
MoFEM::to_non_symm
auto to_non_symm(const FTensor::Tensor2_symmetric< T, DIM > &symm)
Definition Templates.hpp:909
MoFEM::createMFrontInterface
boost::shared_ptr< MFrontInterface > createMFrontInterface(MoFEM::Interface &m_field, ModelHypothesis mh, AssemblyType at)
create mfront interface
Definition MFrontInterface.cpp:2020
MoFEM::ModelHypothesis
ModelHypothesis
Enumeration of model hypotheses supported by MFront interface.
Definition MFrontInterface.hpp:109
MoFEM::AXISYMMETRICAL
@ AXISYMMETRICAL
Axisymmetrical model hypothesis.
Definition MFrontInterface.hpp:112
MoFEM::PLANESTRAIN
@ PLANESTRAIN
Plane strain model hypothesis.
Definition MFrontInterface.hpp:111
MoFEM::TRIDIMENSIONAL
@ TRIDIMENSIONAL
3D model hypothesis.
Definition MFrontInterface.hpp:110
MoFEM::getVoigtVecSymmAxisymm
auto getVoigtVecSymmAxisymm(T &t_mat, const double hoop_term)
Definition Templates.hpp:1919
PlasticOps::opFactoryDomainRhs
MoFEMErrorCode opFactoryDomainRhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
Definition PlasticOps.hpp:503
PlasticOps::opFactoryDomainLhs
MoFEMErrorCode opFactoryDomainLhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
Definition PlasticOps.hpp:550
sdf.r
int r
Definition sdf.py:205
OpGradTimesTensor
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
Definition operators_tests.cpp:48
I
constexpr IntegrationType I
Definition operators_tests.cpp:31
h
double h
Definition photon_diffusion.cpp:71
field_name
constexpr auto field_name
Definition poisson_2d_homogeneous.cpp:13
OpBase
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition radiation.cpp:29
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition EntitiesFieldData.hpp:150
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition EntitiesFieldData.hpp:1806
MoFEM::EntitiesFieldData::EntData::getFTensor1DiffN
auto getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition EntitiesFieldData.hpp:802
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition EntitiesFieldData.hpp:1600
MoFEM::MFrontInterface::CommonData::CommonData
CommonData(MoFEM::Interface &m_field)
Definition MFrontInterface.cpp:2016
MoFEM::OpBaseImpl::nbRows
int nbRows
number of dofs on rows
Definition FormsIntegrators.hpp:368
MoFEM::OpBaseImpl::nbIntegrationPts
int nbIntegrationPts
number of integration points
Definition FormsIntegrators.hpp:370
MoFEM::OpBaseImpl::iNtegrate
virtual MoFEMErrorCode iNtegrate(EntData &, EntData &)
Integrate grad-grad operator.
Definition FormsIntegrators.hpp:391
MoFEM::OpBaseImpl::nbCols
int nbCols
number if dof on column
Definition FormsIntegrators.hpp:369
MoFEM::OpBaseImpl::nbRowBaseFunctions
int nbRowBaseFunctions
number or row base functions
Definition FormsIntegrators.hpp:371
AT
constexpr AssemblyType AT
Definition test_broken_space.cpp:24
is_large_strain
PetscBool is_large_strain
Definition contact.cpp:92
Generated by Doxygen 1.9.8 and hosted at University of Glasgow