v0.15.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("src/libBehaviour.so"),
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 %ld "
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, grad_size);
405 CHKERR mget_tag_data(gradientTag, mPrevGradPtr, stress_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 =
575 "src/libBehaviour." + string(DEFAULT_LIB_EXTENSION);
576 CHKERR PetscOptionsString(
577 param_path.c_str(), "path to the behaviour library", "",
578 default_lib_path.c_str(), char_name, 255, &is_param);
579 if (is_param)
580 lib_path = string(char_name);
581 auto &mgis_bv_ptr = block.second.mGisBehaviour;
582 bool is_finite_strain = false;
583
584 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform)
585 << "Loading behaviour from " << lib_path;
586
587 CHKERR check_lib_finite_strain(lib_path, name, is_finite_strain);
588
589 mgis::behaviour::Hypothesis h;
590 switch (MH) {
591 case TRIDIMENSIONAL:
592 h = mgis::behaviour::Hypothesis::TRIDIMENSIONAL;
593 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform)
594 << "Model hypothesis: TRIDIMENSIONAL";
595 break;
596 case PLANESTRAIN:
597 h = mgis::behaviour::Hypothesis::PLANESTRAIN;
598 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform)
599 << "Model hypothesis: PLANESTRAIN";
600 break;
601 case AXISYMMETRICAL:
602 h = mgis::behaviour::Hypothesis::AXISYMMETRICAL;
603 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform)
604 << "Model hypothesis: AXISYMMETRICAL";
605 break;
606 default:
607 break;
608 }
609
610 if (is_finite_strain) {
611 mgis_bv_ptr = boost::make_shared<Behaviour>(load(op, lib_path, name, h));
612 block.second.isFiniteStrain = true;
613 } else
614 mgis_bv_ptr = boost::make_shared<Behaviour>(load(lib_path, name, h));
615
616 CHKERR block.second.setBlockBehaviourData(set_from_blocks);
617 for (size_t dd = 0; dd < mgis_bv_ptr->mps.size(); ++dd) {
618 double my_param = 0;
619 PetscBool is_set = PETSC_FALSE;
620 string param_cmd = "-param_" + to_string(id) + "_" + to_string(dd);
621 CHKERR PetscOptionsScalar(param_cmd.c_str(), "parameter from cmd", "",
622 my_param, &my_param, &is_set);
623 if (!is_set)
624 continue;
625 setMaterialProperty(block.second.behDataPtr->s0, dd, my_param);
626 setMaterialProperty(block.second.behDataPtr->s1, dd, my_param);
627 }
628
629 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform)
630 << mgis_bv_ptr->behaviour << " behaviour loaded on block "
631 << block.first;
632
633 if (is_finite_strain)
634 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform)
635 << "Finite Strain Kinematics";
636 else
637 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform)
638 << "Small Strain Kinematics";
639
640 if (mgis_bv_ptr->isvs.size()) {
641 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform) << "Internal variables:";
642 for (const auto &is : mgis_bv_ptr->isvs) {
643 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform) << ": " << is.name;
644 }
645 }
646
647 if (mgis_bv_ptr->esvs.size()) {
648 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform) << "External variables:";
649 for (const auto &es : mgis_bv_ptr->esvs) {
650 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform) << ": " << es.name;
651 }
652 }
653
654 auto it = block.second.behDataPtr->s0.material_properties.begin();
655 int nb = 0;
656
657 if (mgis_bv_ptr->mps.size()) {
658 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform) << "Material properties:";
659 for (const auto &mp : mgis_bv_ptr->mps) {
660 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform)
661 << nb++ << " : " << mp.name << " = " << *it++;
662 }
663 }
664
665 if (mgis_bv_ptr->params.size()) {
666 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform) << "Real parameters:";
667 for (const auto &p : mgis_bv_ptr->params) {
668 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform) << nb++ << " : " << p;
669 }
670 }
671
672 if (mgis_bv_ptr->iparams.size()) {
673 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform) << "Integer parameters:";
674 for (const auto &p : mgis_bv_ptr->iparams) {
675 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform) << nb++ << " : " << p;
676 }
677 }
678
679 if (mgis_bv_ptr->usparams.size()) {
680 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform)
681 << "Unsigned short parameters:";
682 for (const auto &p : mgis_bv_ptr->usparams) {
683 MOFEM_LOG("MFrontInterfaceWorld", Sev::inform) << nb++ << " : " << p;
684 }
685 }
686 }
687
688 PetscOptionsEnd();
689
690 auto check_behaviours_kinematics = [&](bool &is_finite_kin) {
691 MoFEMFunctionBeginHot;
692 is_finite_kin =
693 commonDataPtr->setOfBlocksData.begin()->second.isFiniteStrain;
694 for (auto &block : commonDataPtr->setOfBlocksData) {
695 if (block.second.isFiniteStrain != is_finite_kin)
696 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
697 "All used MFront behaviours have to be of same kinematics "
698 "(small or "
699 "large strains)");
700 }
701 MoFEMFunctionReturnHot(0);
702 };
703
704 CHKERR check_behaviours_kinematics(isFiniteKinematics);
705
706 auto get_base = [&]() {
707 Range domain_ents;
708 CHKERR mField.get_moab().get_entities_by_dimension(0, DIM, domain_ents,
709 true);
710 if (domain_ents.empty())
711 CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Empty mesh");
712 const auto type = type_from_handle(domain_ents[0]);
713 switch (type) {
714 case MBQUAD:
715 return DEMKOWICZ_JACOBI_BASE;
716 case MBHEX:
717 return DEMKOWICZ_JACOBI_BASE;
718 case MBTRI:
719 return AINSWORTH_LEGENDRE_BASE;
720 case MBTET:
721 return AINSWORTH_LEGENDRE_BASE;
722 default:
723 CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Element type not handled");
724 }
725 return NOBASE;
726 };
727
728 fieldBase = get_base();
729 MOFEM_LOG("MFrontInterfaceWorld", Sev::verbose)
730 << "MFront post process projection base: "
731 << ApproximationBaseNames[fieldBase];
732
733 MoFEMFunctionReturn(0);
734}
735
736template <int DIM>
737using Tensor4Pack =
738 FTensor::Tensor4<FTensor::PackPtr<double *, 1>, DIM, DIM, DIM, DIM>;
739
740template <int DIM>
741using Tensor2Pack = FTensor::Tensor2<FTensor::PackPtr<double *, 1>, DIM, DIM>;
742
743template <int DIM>
744using Tensor1Pack = FTensor::Tensor1<FTensor::PackPtr<double *, 1>, DIM>;
745
746using Tensor1PackCoords = FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>;
747
748template <int DIM>
749using DdgPack = FTensor::Ddg<FTensor::PackPtr<double *, 1>, DIM, DIM>;
750
751template <typename T> inline size_t get_paraview_size(T &vsize) {
752 return vsize > 1 ? (vsize > 3 ? 9 : 3) : 1;
753};
754
755
756template <typename T1, typename T2>
757inline MoFEMErrorCode get_tensor4_from_voigt(const T1 &K, T2 &D) {
758 MoFEMFunctionBeginHot;
759
760 Number<0> N0;
761 Number<1> N1;
762 Number<2> N2;
763
764 // 3D finite strain
765 if (std::is_same<T2, Tensor4Pack<3>>::value) {
766 D(N0, N0, N0, N0) = K[0];
767 D(N0, N0, N1, N1) = K[1];
768 D(N0, N0, N2, N2) = K[2];
769 D(N0, N0, N0, N1) = K[3];
770 D(N0, N0, N1, N0) = K[4];
771 D(N0, N0, N0, N2) = K[5];
772 D(N0, N0, N2, N0) = K[6];
773 D(N0, N0, N1, N2) = K[7];
774 D(N0, N0, N2, N1) = K[8];
775 D(N1, N1, N0, N0) = K[9];
776 D(N1, N1, N1, N1) = K[10];
777 D(N1, N1, N2, N2) = K[11];
778 D(N1, N1, N0, N1) = K[12];
779 D(N1, N1, N1, N0) = K[13];
780 D(N1, N1, N0, N2) = K[14];
781 D(N1, N1, N2, N0) = K[15];
782 D(N1, N1, N1, N2) = K[16];
783 D(N1, N1, N2, N1) = K[17];
784 D(N2, N2, N0, N0) = K[18];
785 D(N2, N2, N1, N1) = K[19];
786 D(N2, N2, N2, N2) = K[20];
787 D(N2, N2, N0, N1) = K[21];
788 D(N2, N2, N1, N0) = K[22];
789 D(N2, N2, N0, N2) = K[23];
790 D(N2, N2, N2, N0) = K[24];
791 D(N2, N2, N1, N2) = K[25];
792 D(N2, N2, N2, N1) = K[26];
793 D(N0, N1, N0, N0) = K[27];
794 D(N0, N1, N1, N1) = K[28];
795 D(N0, N1, N2, N2) = K[29];
796 D(N0, N1, N0, N1) = K[30];
797 D(N0, N1, N1, N0) = K[31];
798 D(N0, N1, N0, N2) = K[32];
799 D(N0, N1, N2, N0) = K[33];
800 D(N0, N1, N1, N2) = K[34];
801 D(N0, N1, N2, N1) = K[35];
802 D(N1, N0, N0, N0) = K[36];
803 D(N1, N0, N1, N1) = K[37];
804 D(N1, N0, N2, N2) = K[38];
805 D(N1, N0, N0, N1) = K[39];
806 D(N1, N0, N1, N0) = K[40];
807 D(N1, N0, N0, N2) = K[41];
808 D(N1, N0, N2, N0) = K[42];
809 D(N1, N0, N1, N2) = K[43];
810 D(N1, N0, N2, N1) = K[44];
811 D(N0, N2, N0, N0) = K[45];
812 D(N0, N2, N1, N1) = K[46];
813 D(N0, N2, N2, N2) = K[47];
814 D(N0, N2, N0, N1) = K[48];
815 D(N0, N2, N1, N0) = K[49];
816 D(N0, N2, N0, N2) = K[50];
817 D(N0, N2, N2, N0) = K[51];
818 D(N0, N2, N1, N2) = K[52];
819 D(N0, N2, N2, N1) = K[53];
820 D(N2, N0, N0, N0) = K[54];
821 D(N2, N0, N1, N1) = K[55];
822 D(N2, N0, N2, N2) = K[56];
823 D(N2, N0, N0, N1) = K[57];
824 D(N2, N0, N1, N0) = K[58];
825 D(N2, N0, N0, N2) = K[59];
826 D(N2, N0, N2, N0) = K[60];
827 D(N2, N0, N1, N2) = K[61];
828 D(N2, N0, N2, N1) = K[62];
829 D(N1, N2, N0, N0) = K[63];
830 D(N1, N2, N1, N1) = K[64];
831 D(N1, N2, N2, N2) = K[65];
832 D(N1, N2, N0, N1) = K[66];
833 D(N1, N2, N1, N0) = K[67];
834 D(N1, N2, N0, N2) = K[68];
835 D(N1, N2, N2, N0) = K[69];
836 D(N1, N2, N1, N2) = K[70];
837 D(N1, N2, N2, N1) = K[71];
838 D(N2, N1, N0, N0) = K[72];
839 D(N2, N1, N1, N1) = K[73];
840 D(N2, N1, N2, N2) = K[74];
841 D(N2, N1, N0, N1) = K[75];
842 D(N2, N1, N1, N0) = K[76];
843 D(N2, N1, N0, N2) = K[77];
844 D(N2, N1, N2, N0) = K[78];
845 D(N2, N1, N1, N2) = K[79];
846 D(N2, N1, N2, N1) = K[80];
847 }
848
849 // plane strain, finite strain
850 if (std::is_same<T2, Tensor4Pack<2>>::value) {
851 D(N0, N0, N0, N0) = K[0];
852 D(N0, N0, N1, N1) = K[1];
853 // K[2] is not used
854 D(N0, N0, N0, N1) = K[3];
855 D(N0, N0, N1, N0) = K[4];
856 D(N1, N1, N0, N0) = K[5];
857 D(N1, N1, N1, N1) = K[6];
858 // K[7] is not used
859 D(N1, N1, N0, N1) = K[8];
860 D(N1, N1, N1, N0) = K[9];
861 // K[10], K[11], K[12], K[13], K[14] are not used
862 D(N0, N1, N0, N0) = K[15];
863 D(N0, N1, N1, N1) = K[16];
864 // K[17] is not used
865 D(N0, N1, N0, N1) = K[18];
866 D(N0, N1, N1, N0) = K[19];
867 D(N1, N0, N0, N0) = K[20];
868 D(N1, N0, N1, N1) = K[21];
869 // K[22] is not used
870 D(N1, N0, N0, N1) = K[23];
871 D(N1, N0, N1, N0) = K[24];
872 }
873
874 // 3D small strain
875 if (std::is_same<T2, DdgPack<3>>::value) {
876 D(N0, N0, N0, N0) = K[0];
877 D(N0, N0, N1, N1) = K[1];
878 D(N0, N0, N2, N2) = K[2];
879
880 D(N0, N0, N0, N1) = inv_sqr2 * K[3];
881 D(N0, N0, N0, N2) = inv_sqr2 * K[4];
882 D(N0, N0, N1, N2) = inv_sqr2 * K[5];
883
884 D(N1, N1, N0, N0) = K[6];
885 D(N1, N1, N1, N1) = K[7];
886 D(N1, N1, N2, N2) = K[8];
887
888 D(N1, N1, N0, N1) = inv_sqr2 * K[9];
889 D(N1, N1, N0, N2) = inv_sqr2 * K[10];
890 D(N1, N1, N1, N2) = inv_sqr2 * K[11];
891
892 D(N2, N2, N0, N0) = K[12];
893 D(N2, N2, N1, N1) = K[13];
894 D(N2, N2, N2, N2) = K[14];
895
896 D(N2, N2, N0, N1) = inv_sqr2 * K[15];
897 D(N2, N2, N0, N2) = inv_sqr2 * K[16];
898 D(N2, N2, N1, N2) = inv_sqr2 * K[17];
899
900 D(N0, N1, N0, N0) = inv_sqr2 * K[18];
901 D(N0, N1, N1, N1) = inv_sqr2 * K[19];
902 D(N0, N1, N2, N2) = inv_sqr2 * K[20];
903
904 D(N0, N1, N0, N1) = 0.5 * K[21];
905 D(N0, N1, N0, N2) = 0.5 * K[22];
906 D(N0, N1, N1, N2) = 0.5 * K[23];
907
908 D(N0, N2, N0, N0) = inv_sqr2 * K[24];
909 D(N0, N2, N1, N1) = inv_sqr2 * K[25];
910 D(N0, N2, N2, N2) = inv_sqr2 * K[26];
911
912 D(N0, N2, N0, N1) = 0.5 * K[27];
913 D(N0, N2, N0, N2) = 0.5 * K[28];
914 D(N0, N2, N1, N2) = 0.5 * K[29];
915
916 D(N1, N2, N0, N0) = inv_sqr2 * K[30];
917 D(N1, N2, N1, N1) = inv_sqr2 * K[31];
918 D(N1, N2, N2, N2) = inv_sqr2 * K[32];
919
920 D(N1, N2, N0, N1) = 0.5 * K[33];
921 D(N1, N2, N0, N2) = 0.5 * K[34];
922 D(N1, N2, N1, N2) = 0.5 * K[35];
923 }
924
925 // plane strain, small strain
926 if (std::is_same<T2, DdgPack<2>>::value) {
927
928 D(N0, N0, N0, N0) = K[0];
929 D(N0, N0, N1, N1) = K[1];
930 // K[2] is not used
931
932 D(N0, N0, N0, N1) = inv_sqr2 * K[3];
933
934 D(N1, N1, N0, N0) = K[4];
935 D(N1, N1, N1, N1) = K[5];
936 // K[6] is not used
937
938 D(N1, N1, N0, N1) = inv_sqr2 * K[7];
939
940 // K[8], K[9], K[10], K[11] are not used
941
942 D(N0, N1, N0, N0) = inv_sqr2 * K[12];
943 D(N0, N1, N1, N1) = inv_sqr2 * K[13];
944
945 // K[14] is not used
946
947 D(N0, N1, N0, N1) = 0.5 * K[15];
948 }
949
950 MoFEMFunctionReturnHot(0);
951};
952
953template <bool IS_LARGE_STRAIN, typename T1, typename T2>
954inline MoFEMErrorCode get_full_tensor4_from_voigt(const T1 &K, T2 &D) {
955 MoFEMFunctionBeginHot;
956
957 Number<0> N0;
958 Number<1> N1;
959 Number<2> N2;
960
961 // 2D finite strain, full tensor
962 if constexpr (IS_LARGE_STRAIN) {
963 D(N0, N0, N0, N0) = K[0];
964 D(N0, N0, N1, N1) = K[1];
965 D(N0, N0, N2, N2) = K[2];
966 D(N0, N0, N0, N1) = K[3];
967 D(N0, N0, N1, N0) = K[4];
968 D(N1, N1, N0, N0) = K[5];
969 D(N1, N1, N1, N1) = K[6];
970 D(N1, N1, N2, N2) = K[7];
971 D(N1, N1, N0, N1) = K[8];
972 D(N1, N1, N1, N0) = K[9];
973 D(N2, N2, N0, N0) = K[10];
974 D(N2, N2, N1, N1) = K[11];
975 D(N2, N2, N2, N2) = K[12];
976 D(N2, N2, N0, N1) = K[13];
977 D(N2, N2, N1, N0) = K[14];
978 D(N0, N1, N0, N0) = K[15];
979 D(N0, N1, N1, N1) = K[16];
980 D(N0, N1, N2, N2) = K[17];
981 D(N0, N1, N0, N1) = K[18];
982 D(N0, N1, N1, N0) = K[19];
983 D(N1, N0, N0, N0) = K[20];
984 D(N1, N0, N1, N1) = K[21];
985 D(N1, N0, N2, N2) = K[22];
986 D(N1, N0, N0, N1) = K[23];
987 D(N1, N0, N1, N0) = K[24];
988 } else { // 2D small strain, full tensor
989
990 D(N0, N0, N0, N0) = K[0];
991 D(N0, N0, N1, N1) = K[1];
992 D(N0, N0, N2, N2) = K[2];
993
994 D(N0, N0, N0, N1) = inv_sqr2 * K[3];
995
996 D(N1, N1, N0, N0) = K[4];
997 D(N1, N1, N1, N1) = K[5];
998 D(N1, N1, N2, N2) = K[6];
999
1000 D(N1, N1, N0, N1) = inv_sqr2 * K[7];
1001
1002 D(N2, N2, N0, N0) = K[8];
1003 D(N2, N2, N1, N1) = K[9];
1004 D(N2, N2, N2, N2) = K[10];
1005
1006 D(N2, N2, N0, N1) = inv_sqr2 * K[11];
1007 D(N2, N2, N1, N0) = D(N2, N2, N0, N1);
1008
1009 D(N0, N1, N0, N0) = inv_sqr2 * K[12];
1010 D(N0, N1, N1, N1) = inv_sqr2 * K[13];
1011
1012 D(N0, N1, N2, N2) = inv_sqr2 * K[14];
1013 D(N1, N0, N2, N2) = D(N0, N1, N2, N2);
1014
1015 D(N0, N1, N0, N1) = 0.5 * K[15];
1016 }
1017
1018 MoFEMFunctionReturnHot(0);
1019};
1020
1021template <typename T> T get_tangent_tensor(MatrixDouble &mat);
1022
1023template <>
1024Tensor4Pack<3> get_tangent_tensor<Tensor4Pack<3>>(MatrixDouble &mat) {
1025 return getFTensor4FromMat<3, 3, 3, 3>(mat);
1026}
1027
1028template <>
1029Tensor4Pack<2> get_tangent_tensor<Tensor4Pack<2>>(MatrixDouble &mat) {
1030 return getFTensor4FromMat<2, 2, 2, 2>(mat);
1031}
1032
1033template <> DdgPack<3> get_tangent_tensor<DdgPack<3>>(MatrixDouble &mat) {
1034 return getFTensor4DdgFromMat<3, 3>(mat);
1035}
1036
1037template <> DdgPack<2> get_tangent_tensor<DdgPack<2>>(MatrixDouble &mat) {
1038 return getFTensor4DdgFromMat<2, 2>(mat);
1039}
1040
1041template <bool IS_LARGE_STRAIN, ModelHypothesis MH>
1042MoFEMErrorCode
1043mgis_integration(size_t gg, Tensor2Pack<MFrontEleType<MH>::SPACE_DIM> &t_grad,
1044 Tensor1Pack<MFrontEleType<MH>::SPACE_DIM> &t_disp,
1045 Tensor1PackCoords &t_coords,
1046 MFrontInterface::CommonData &common_data,
1047 MFrontInterface::CommonData::BlockData &block_data) {
1048 MoFEMFunctionBegin;
1049
1050 static constexpr int DIM = MFrontEleType<MH>::SPACE_DIM;
1051
1052 int check_integration;
1053 MatrixDouble &mat_int = *common_data.internalVariablePtr;
1054 MatrixDouble &mat_grad0 = *common_data.mPrevGradPtr;
1055 MatrixDouble &mat_stress0 = *common_data.mPrevStressPtr;
1056
1057 int &size_of_vars = block_data.sizeIntVar;
1058 int &size_of_grad = block_data.sizeGradVar;
1059 int &size_of_stress = block_data.sizeStressVar;
1060
1061 auto &mgis_bv = *block_data.mGisBehaviour;
1062
1063 FTensor::Index<'i', DIM> i;
1064 FTensor::Index<'j', DIM> j;
1065
1066 if constexpr (IS_LARGE_STRAIN) {
1067 FTensor::Tensor2<double, DIM, DIM> t_strain;
1068 t_strain(i, j) = t_grad(i, j) + kronecker_delta(i, j);
1069 if constexpr (MH == AXISYMMETRICAL) {
1070 setGradient(block_data.behDataPtr->s1, 0, size_of_grad,
1071 &*getVoigtVecAxisymm(t_strain, t_disp(0) / t_coords(0)).data());
1072 } else {
1073 setGradient(block_data.behDataPtr->s1, 0, size_of_grad,
1074 &*getVoigtVec<DIM>(t_strain).data());
1075 }
1076 } else {
1077 FTensor::Tensor2_symmetric<double, DIM> t_strain;
1078 t_strain(i, j) = (t_grad(i, j) || t_grad(j, i)) / 2;
1079 if constexpr (MH == AXISYMMETRICAL) {
1080 setGradient(
1081 block_data.behDataPtr->s1, 0, size_of_grad,
1082 &*getVoigtVecSymmAxisymm(t_strain, t_disp(0) / t_coords(0)).data());
1083 } else {
1084 setGradient(block_data.behDataPtr->s1, 0, size_of_grad,
1085 &*getVoigtVecSymm<DIM>(t_strain).data());
1086 }
1087 }
1088
1089 auto grad0_vec =
1090 getVectorAdaptor(&mat_grad0.data()[gg * size_of_grad], size_of_grad);
1091 setGradient(block_data.behDataPtr->s0, 0, size_of_grad, &*grad0_vec.begin());
1092
1093 auto stress0_vec = getVectorAdaptor(&mat_stress0.data()[gg * size_of_stress],
1094 size_of_stress);
1095 setThermodynamicForce(block_data.behDataPtr->s0, 0, size_of_stress,
1096 &*stress0_vec.begin());
1097
1098 if (size_of_vars) {
1099 auto internal_var =
1100 getVectorAdaptor(&mat_int.data()[gg * size_of_vars], size_of_vars);
1101 setInternalStateVariable(block_data.behDataPtr->s0, 0, size_of_vars,
1102 &*internal_var.begin());
1103 }
1104
1105 check_integration = mgis::behaviour::integrate(block_data.bView, mgis_bv);
1106 switch (check_integration) {
1107 case -1:
1108 MOFEM_LOG("WORLD", Sev::error) << "Mfront integration failed";
1109 break;
1110 case 0:
1111 MOFEM_LOG("WORLD", Sev::warning)
1112 << "Mfront integration succeeded but results are unreliable";
1113 break;
1114 case 1:
1115 default:
1116 break;
1117 }
1118
1119 MoFEMFunctionReturn(0);
1120}
1121
1122template <bool UPDATE, bool IS_LARGE_STRAIN, ModelHypothesis MH>
1123struct OpStressTmp : public MFrontEleType<MH>::DomainEleOp {
1124 static constexpr int DIM = MFrontEleType<MH>::SPACE_DIM;
1125 using DomainEleOp = typename MFrontEleType<MH>::DomainEleOp;
1126
1127 OpStressTmp(const std::string field_name,
1128 boost::shared_ptr<MFrontInterface::CommonData> common_data_ptr,
1129 boost::shared_ptr<FEMethod> monitor_ptr)
1130 : DomainEleOp(field_name, DomainEleOp::OPROW),
1131 commonDataPtr(common_data_ptr), monitorPtr(monitor_ptr) {
1132 std::fill(&DomainEleOp::doEntities[MBEDGE],
1133 &DomainEleOp::doEntities[MBMAXTYPE], false);
1134 }
1135 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
1136
1137private:
1138 boost::shared_ptr<MFrontInterface::CommonData> commonDataPtr;
1139 boost::shared_ptr<FEMethod> monitorPtr;
1140};
1141
1142template <ModelHypothesis MH>
1143using OpUpdateVariablesFiniteStrains = OpStressTmp<true, true, MH>;
1144
1145template <ModelHypothesis MH>
1146using OpUpdateVariablesSmallStrains = OpStressTmp<true, false, MH>;
1147
1148template <ModelHypothesis MH>
1149using OpStressFiniteStrains = OpStressTmp<false, true, MH>;
1150
1151template <ModelHypothesis MH>
1152using OpStressSmallStrains = OpStressTmp<false, false, MH>;
1153
1154template <bool UPDATE, bool IS_LARGE_STRAIN, ModelHypothesis MH>
1155MoFEMErrorCode OpStressTmp<UPDATE, IS_LARGE_STRAIN, MH>::doWork(int side,
1156 EntityType type,
1157 EntData &data) {
1158 MoFEMFunctionBegin;
1159
1160 Index<'i', 3> i;
1161 Index<'j', 3> j;
1162
1163 Index<'I', 2> I;
1164 Index<'J', 2> J;
1165
1166 const size_t nb_gauss_pts = commonDataPtr->mGradPtr->size2();
1167 auto fe_ent = DomainEleOp::getNumeredEntFiniteElementPtr()->getEnt();
1168 auto id = commonDataPtr->blocksIDmap.at(fe_ent);
1169 auto &dAta = commonDataPtr->setOfBlocksData.at(id);
1170
1171 if (monitorPtr == nullptr)
1172 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_FOUND,
1173 "Time Monitor (FEMethod) has not been set for MFrontInterfaceImpl. "
1174 "Make sure to call setMonitor before calling "
1175 "opFactoryDomainRhs and opFactoryDomainLhs");
1176
1177 dAta.setTag(DataTags::RHS);
1178 dAta.behDataPtr->dt = monitorPtr->ts_dt;
1179 dAta.bView.dt = monitorPtr->ts_dt;
1180
1181 CHKERR commonDataPtr->getInternalVar(fe_ent, nb_gauss_pts, dAta.sizeIntVar,
1182 dAta.sizeGradVar, dAta.sizeStressVar,
1183 IS_LARGE_STRAIN);
1184
1185 MatrixDouble &mat_int = *commonDataPtr->internalVariablePtr;
1186 MatrixDouble &mat_grad0 = *commonDataPtr->mPrevGradPtr;
1187 MatrixDouble &mat_stress0 = *commonDataPtr->mPrevStressPtr;
1188
1189 auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->mGradPtr));
1190 auto t_disp = getFTensor1FromMat<DIM>(*(commonDataPtr->mDispPtr));
1191 auto t_coords = DomainEleOp::getFTensor1CoordsAtGaussPts();
1192
1193 commonDataPtr->mStressPtr->resize(DIM * DIM, nb_gauss_pts);
1194 commonDataPtr->mStressPtr->clear();
1195 auto t_stress = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->mStressPtr));
1196
1197 commonDataPtr->mFullStressPtr->resize(3 * 3, nb_gauss_pts);
1198 commonDataPtr->mFullStressPtr->clear();
1199 auto t_full_stress =
1200 getFTensor2FromMat<3, 3>(*(commonDataPtr->mFullStressPtr));
1201
1202 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1203
1204 CHKERR mgis_integration<IS_LARGE_STRAIN, MH>(gg, t_grad, t_disp, t_coords,
1205 *commonDataPtr, dAta);
1206
1207 if constexpr (DIM == 3) {
1208 Tensor2<double, 3, 3> t_force;
1209 if constexpr (IS_LARGE_STRAIN) {
1210 t_force = Tensor2<double, 3, 3>(
1211 VOIGT_VEC_3D(getThermodynamicForce(dAta.behDataPtr->s1, 0)));
1212 } else {
1213 t_force = to_non_symm(Tensor2_symmetric<double, 3>(
1214 VOIGT_VEC_SYMM_3D(getThermodynamicForce(dAta.behDataPtr->s1, 0))));
1215 }
1216 t_stress(i, j) = t_force(i, j);
1217 } else if constexpr (DIM == 2) {
1218 Tensor2<double, 2, 2> t_force;
1219 Tensor2<double, 3, 3> t_full_force;
1220 if constexpr (IS_LARGE_STRAIN) {
1221 t_force = Tensor2<double, 2, 2>(
1222 VOIGT_VEC_2D(getThermodynamicForce(dAta.behDataPtr->s1, 0)));
1223 t_full_force = Tensor2<double, 3, 3>(
1224 VOIGT_VEC_2D_FULL(getThermodynamicForce(dAta.behDataPtr->s1, 0)));
1225 } else {
1226 t_force = to_non_symm(Tensor2_symmetric<double, 2>(
1227 VOIGT_VEC_SYMM_2D(getThermodynamicForce(dAta.behDataPtr->s1, 0))));
1228 t_full_force =
1229 to_non_symm(Tensor2_symmetric<double, 3>(VOIGT_VEC_SYMM_2D_FULL(
1230 getThermodynamicForce(dAta.behDataPtr->s1, 0))));
1231 }
1232 t_stress(I, J) = t_force(I, J);
1233 t_full_stress(i, j) = t_full_force(i, j);
1234 }
1235
1236 if constexpr (UPDATE) {
1237 for (int dd = 0; dd != dAta.sizeIntVar; ++dd) {
1238 mat_int(gg, dd) = *getInternalStateVariable(dAta.behDataPtr->s1, dd);
1239 }
1240 for (int dd = 0; dd != dAta.sizeGradVar; ++dd) {
1241 mat_grad0(gg, dd) = *getGradient(dAta.behDataPtr->s1, dd);
1242 }
1243 for (int dd = 0; dd != dAta.sizeStressVar; ++dd) {
1244 mat_stress0(gg, dd) = *getThermodynamicForce(dAta.behDataPtr->s1, dd);
1245 }
1246 }
1247
1248 ++t_stress;
1249 ++t_full_stress;
1250 ++t_grad;
1251 ++t_disp;
1252 ++t_coords;
1253 }
1254
1255 if constexpr (UPDATE) {
1256 CHKERR commonDataPtr->setInternalVar(fe_ent);
1257 }
1258
1259 MoFEMFunctionReturn(0);
1260}
1261
1262template <typename T, ModelHypothesis MH>
1263struct OpTangent : public MFrontEleType<MH>::DomainEleOp {
1264 static constexpr int DIM = MFrontEleType<MH>::SPACE_DIM;
1265 using DomainEleOp = typename MFrontEleType<MH>::DomainEleOp;
1266
1267 OpTangent(const std::string field_name,
1268 boost::shared_ptr<MFrontInterface::CommonData> common_data_ptr,
1269 boost::shared_ptr<FEMethod> monitor_ptr)
1270 : DomainEleOp(field_name, DomainEleOp::OPROW),
1271 commonDataPtr(common_data_ptr), monitorPtr(monitor_ptr) {
1272 std::fill(&DomainEleOp::doEntities[MBEDGE],
1273 &DomainEleOp::doEntities[MBMAXTYPE], false);
1274 }
1275 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
1276
1277private:
1278 boost::shared_ptr<MFrontInterface::CommonData> commonDataPtr;
1279 boost::shared_ptr<FEMethod> monitorPtr;
1280};
1281
1282template <ModelHypothesis MH>
1283using OpTangentFiniteStrains =
1284 struct OpTangent<Tensor4Pack<MFrontEleType<MH>::SPACE_DIM>, MH>;
1285
1286template <ModelHypothesis MH>
1287using OpTangentSmallStrains =
1288 struct OpTangent<DdgPack<MFrontEleType<MH>::SPACE_DIM>, MH>;
1289
1290template <typename T, ModelHypothesis MH>
1291MoFEMErrorCode OpTangent<T, MH>::doWork(int side, EntityType type,
1292 EntData &data) {
1293 MoFEMFunctionBegin;
1294
1295 const size_t nb_gauss_pts = commonDataPtr->mGradPtr->size2();
1296 auto fe_ent = DomainEleOp::getNumeredEntFiniteElementPtr()->getEnt();
1297 auto id = commonDataPtr->blocksIDmap.at(fe_ent);
1298 auto &dAta = commonDataPtr->setOfBlocksData.at(id);
1299
1300 if (monitorPtr == nullptr)
1301 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_FOUND,
1302 "Time Monitor (FEMethod) has not been set for MFrontInterfaceImpl. "
1303 "Make sure to call setMonitor(monitor_ptr) before calling "
1304 "setOperators().");
1305
1306 dAta.setTag(DataTags::LHS);
1307 dAta.behDataPtr->dt = monitorPtr->ts_dt;
1308 dAta.bView.dt = monitorPtr->ts_dt;
1309
1310 constexpr bool IS_LARGE_STRAIN = std::is_same<T, Tensor4Pack<3>>::value ||
1311 std::is_same<T, Tensor4Pack<2>>::value;
1312
1313 CHKERR commonDataPtr->getInternalVar(fe_ent, nb_gauss_pts, dAta.sizeIntVar,
1314 dAta.sizeGradVar, dAta.sizeStressVar,
1315 IS_LARGE_STRAIN);
1316
1317 MatrixDouble &S_E = *(commonDataPtr->materialTangentPtr);
1318 MatrixDouble &F_E = *(commonDataPtr->mFullTangentPtr);
1319
1320 size_t tens_size = 36;
1321
1322 if constexpr (DIM == 2) {
1323 // plane strain
1324 if constexpr (IS_LARGE_STRAIN)
1325 tens_size = 16;
1326 else
1327 tens_size = 9;
1328 } else {
1329 if constexpr (IS_LARGE_STRAIN) // for finite strains
1330 tens_size = 81;
1331 }
1332
1333 S_E.resize(tens_size, nb_gauss_pts, false);
1334 auto D1 = get_tangent_tensor<T>(S_E);
1335
1336 size_t full_tens_size = 81;
1337 F_E.resize(full_tens_size, nb_gauss_pts, false);
1338 F_E.clear();
1339
1340 auto D2 = get_tangent_tensor<Tensor4Pack<3>>(F_E);
1341
1342 auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->mGradPtr));
1343 auto t_disp = getFTensor1FromMat<DIM>(*(commonDataPtr->mDispPtr));
1344 auto t_coords = DomainEleOp::getFTensor1CoordsAtGaussPts();
1345
1346 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1347
1348 CHKERR mgis_integration<IS_LARGE_STRAIN, MH>(gg, t_grad, t_disp, t_coords,
1349 *commonDataPtr, dAta);
1350
1351 CHKERR get_tensor4_from_voigt(&*dAta.behDataPtr->K.begin(), D1);
1352 CHKERR get_full_tensor4_from_voigt<IS_LARGE_STRAIN>(
1353 &*dAta.behDataPtr->K.begin(), D2);
1354
1355 ++D1;
1356 ++D2;
1357 ++t_grad;
1358 ++t_disp;
1359 ++t_coords;
1360 }
1361
1362 MoFEMFunctionReturn(0);
1363}
1364
1365template <AssemblyType AT>
1366struct OpAxisymmetricRhs
1367 : public FormsIntegrators<
1368 MFrontEleType<AXISYMMETRICAL>::DomainEleOp>::Assembly<AT>::OpBase {
1369 using DomainEleOp = typename MFrontEleType<AXISYMMETRICAL>::DomainEleOp;
1370 using OpBase = typename FormsIntegrators<DomainEleOp>::Assembly<AT>::OpBase;
1371
1372 OpAxisymmetricRhs(
1373 const std::string field_name,
1374 boost::shared_ptr<MFrontInterface::CommonData> common_data_ptr)
1375 : OpBase(field_name, field_name, DomainEleOp::OPROW),
1376 commonDataPtr(common_data_ptr) {};
1377
1378private:
1379 boost::shared_ptr<MFrontInterface::CommonData> commonDataPtr;
1380
1381 MoFEMErrorCode iNtegrate(EntData &row_data);
1382};
1383
1384template <AssemblyType AT>
1385MoFEMErrorCode OpAxisymmetricRhs<AT>::iNtegrate(EntData &row_data) {
1386 MoFEMFunctionBegin;
1387
1388 // get element volume
1389 const double vol = OpBase::getMeasure();
1390 // get integration weights
1391 auto t_w = OpBase::getFTensor0IntegrationWeight();
1392 // get coordinate at integration points
1393
1394 auto t_full_stress =
1395 getFTensor2FromMat<3, 3>(*(commonDataPtr->mFullStressPtr));
1396
1397 // loop over integration points
1398 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
1399
1400 auto t_nf = OpBase::template getNf<2>();
1401
1402 FTensor::Tensor0<double *> t_base(&row_data.getN()(gg, 0));
1403
1404 // take into account Jacobian
1405 const double alpha = t_w * vol * 2. * M_PI;
1406 // loop over rows base functions
1407 for (int rr = 0; rr != OpBase::nbRows / 2; ++rr) {
1408
1409 t_nf(0) += alpha * t_full_stress(2, 2) * t_base;
1410
1411 ++t_base;
1412 ++t_nf;
1413 }
1414
1415 ++t_full_stress;
1416 ++t_w; // move to another integration weight
1417 }
1418
1419 MoFEMFunctionReturn(0);
1420}
1421
1422template <AssemblyType AT>
1423struct OpAxisymmetricLhs
1424 : public FormsIntegrators<
1425 MFrontEleType<AXISYMMETRICAL>::DomainEleOp>::Assembly<AT>::OpBase {
1426 using DomainEleOp = typename MFrontEleType<AXISYMMETRICAL>::DomainEleOp;
1427 using OpBase = typename FormsIntegrators<DomainEleOp>::Assembly<AT>::OpBase;
1428
1429 OpAxisymmetricLhs(
1430 const std::string field_name,
1431 boost::shared_ptr<MFrontInterface::CommonData> common_data_ptr)
1432 : OpBase(field_name, field_name, DomainEleOp::OPROWCOL),
1433 commonDataPtr(common_data_ptr) {};
1434
1435private:
1436 boost::shared_ptr<MFrontInterface::CommonData> commonDataPtr;
1437
1438 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
1439};
1440
1441template <AssemblyType AT>
1442MoFEMErrorCode OpAxisymmetricLhs<AT>::iNtegrate(EntData &row_data,
1443 EntData &col_data) {
1444 MoFEMFunctionBegin;
1445
1446 // get element volume
1447 const double vol = OpBase::getMeasure();
1448 // get integration weights
1449 auto t_w = OpBase::getFTensor0IntegrationWeight();
1450 // get base function gradient on rows
1451 auto t_row_base = row_data.getFTensor0N();
1452 // get derivatives of base functions on rows
1453 auto t_row_diff_base = row_data.getFTensor1DiffN<2>();
1454 // get coordinate at integration points
1455 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
1456
1457 Number<0> N0;
1458 Number<1> N1;
1459 Number<2> N2;
1460
1461 auto t_D =
1462 getFTensor4FromMat<3, 3, 3, 3, 1>(*(commonDataPtr->mFullTangentPtr));
1463
1464 // loop over integration points
1465 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
1466
1467 // Cylinder radius
1468 const double r_cylinder = t_coords(0);
1469
1470 // take into account Jacobean
1471 const double alpha = t_w * vol * 2. * M_PI;
1472
1473 // loop over rows base functions
1474 int rr = 0;
1475 for (; rr != OpBase::nbRows / 2; ++rr) {
1476
1477 // get sub matrix for the row
1478 auto t_m = OpBase::template getLocMat<2>(2 * rr);
1479
1480 // get derivatives of base functions for columns
1481 auto t_col_diff_base = col_data.getFTensor1DiffN<2>(gg, 0);
1482 // get base functions for columns
1483 auto t_col_base = col_data.getFTensor0N(gg, 0);
1484
1485 // loop over columns
1486 for (int cc = 0; cc != OpBase::nbCols / 2; cc++) {
1487
1488 t_m(0, 0) +=
1489 alpha * t_D(N0, N0, N2, N2) * t_col_base * t_row_diff_base(0);
1490
1491 t_m(1, 0) +=
1492 alpha * t_D(N1, N1, N2, N2) * t_col_base * t_row_diff_base(1);
1493
1494 t_m(0, 0) +=
1495 alpha * t_D(N2, N2, N0, N0) * t_col_diff_base(0) * t_row_base;
1496
1497 t_m(0, 1) +=
1498 alpha * t_D(N2, N2, N1, N1) * t_col_diff_base(1) * t_row_base;
1499
1500 t_m(0, 0) += alpha * (t_D(N2, N2, N2, N2) / r_cylinder) * t_col_base *
1501 t_row_base;
1502
1503 t_m(0, 0) +=
1504 alpha * t_D(N2, N2, N0, N1) * t_col_diff_base(1) * t_row_base;
1505
1506 t_m(0, 1) +=
1507 alpha * t_D(N2, N2, N1, N0) * t_col_diff_base(0) * t_row_base;
1508
1509 t_m(0, 0) +=
1510 alpha * t_D(N0, N1, N2, N2) * t_col_base * t_row_diff_base(1);
1511
1512 t_m(1, 0) +=
1513 alpha * t_D(N1, N0, N2, N2) * t_col_base * t_row_diff_base(0);
1514
1515 ++t_col_base;
1516 ++t_col_diff_base;
1517 ++t_m;
1518 }
1519
1520 ++t_row_base;
1521 ++t_row_diff_base;
1522 }
1523
1524 for (; rr < OpBase::nbRowBaseFunctions; ++rr) {
1525 ++t_row_base;
1526 ++t_row_diff_base;
1527 }
1528
1529 ++t_coords;
1530 ++t_w; // move to another integration weight
1531 ++t_D;
1532 }
1533
1534 MoFEMFunctionReturn(0);
1535}
1536
1537template <ModelHypothesis MH, AssemblyType AT>
1538MoFEMErrorCode MFrontInterfaceImpl<MH, AT>::opFactoryDomainRhs(
1539 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1540 std::string field_name) {
1541 MoFEMFunctionBegin;
1542
1543 auto jacobian = [&](const double r, const double, const double) {
1544 if (MH == AXISYMMETRICAL)
1545 return 2. * M_PI * r;
1546 else
1547 return 1.;
1548 };
1549
1550 auto add_domain_ops_rhs = [&](auto &pipeline) {
1551 if (isFiniteKinematics)
1552 pipeline.push_back(
1553 new OpStressFiniteStrains<MH>(field_name, commonDataPtr, monitorPtr));
1554 else
1555 pipeline.push_back(
1556 new OpStressSmallStrains<MH>(field_name, commonDataPtr, monitorPtr));
1557
1558 pipeline.push_back(
1559 new OpInternalForce(field_name, commonDataPtr->mStressPtr, jacobian));
1560
1561 if (MH == AXISYMMETRICAL)
1562 pipeline.push_back(new OpAxisymmetricRhs<AT>(field_name, commonDataPtr));
1563 };
1564
1565 auto add_domain_base_ops = [&](auto &pipeline) {
1566 pipeline.push_back(new OpCalculateVectorFieldValues<DIM>(
1567 field_name, commonDataPtr->mDispPtr));
1568 pipeline.push_back(new OpCalculateVectorFieldGradient<DIM, DIM>(
1569 field_name, commonDataPtr->mGradPtr));
1570 };
1571
1572 add_domain_base_ops(pip);
1573 add_domain_ops_rhs(pip);
1574
1575 MoFEMFunctionReturn(0);
1576}
1577
1578template <ModelHypothesis MH, AssemblyType AT>
1579MoFEMErrorCode MFrontInterfaceImpl<MH, AT>::opFactoryDomainLhs(
1580 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1581 std::string field_name) {
1582 MoFEMFunctionBegin;
1583
1584 auto jacobian = [&](const double r, const double, const double) {
1585 if (MH == AXISYMMETRICAL)
1586 return 2. * M_PI * r;
1587 else
1588 return 1.;
1589 };
1590
1591 auto add_domain_ops_lhs = [&](auto &pipeline) {
1592 if (isFiniteKinematics) {
1593 pipeline.push_back(new OpTangentFiniteStrains<MH>(
1594 field_name, commonDataPtr, monitorPtr));
1595 pipeline.push_back(new OpAssembleLhsFiniteStrains(
1596 field_name, field_name, commonDataPtr->materialTangentPtr, jacobian));
1597 } else {
1598 pipeline.push_back(
1599 new OpTangentSmallStrains<MH>(field_name, commonDataPtr, monitorPtr));
1600 pipeline.push_back(new OpAssembleLhsSmallStrains(
1601 field_name, field_name, commonDataPtr->materialTangentPtr, nullptr,
1602 jacobian));
1603 }
1604 if (MH == AXISYMMETRICAL)
1605 pipeline.push_back(new OpAxisymmetricLhs<AT>(field_name, commonDataPtr));
1606 };
1607
1608 auto add_domain_base_ops = [&](auto &pipeline) {
1609 pipeline.push_back(new OpCalculateVectorFieldValues<DIM>(
1610 field_name, commonDataPtr->mDispPtr));
1611 pipeline.push_back(new OpCalculateVectorFieldGradient<DIM, DIM>(
1612 field_name, commonDataPtr->mGradPtr));
1613 };
1614
1615 add_domain_base_ops(pip);
1616 add_domain_ops_lhs(pip);
1617
1618 MoFEMFunctionReturn(0);
1619}
1620
1621template <ModelHypothesis MH>
1622struct OpSaveGaussPts : public MFrontEleType<MH>::DomainEleOp {
1623 static constexpr int DIM = MFrontEleType<MH>::SPACE_DIM;
1624 using DomainEleOp = typename MFrontEleType<MH>::DomainEleOp;
1625
1626 OpSaveGaussPts(const std::string field_name, moab::Interface &moab_mesh,
1627 boost::shared_ptr<MFrontInterface::CommonData> common_data_ptr)
1628 : DomainEleOp(field_name, DomainEleOp::OPROW), internalVarMesh(moab_mesh),
1629 commonDataPtr(common_data_ptr), fieldName(field_name) {
1630 std::fill(&DomainEleOp::doEntities[MBEDGE],
1631 &DomainEleOp::doEntities[MBMAXTYPE], false);
1632 }
1633
1634 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1635 MoFEMFunctionBegin;
1636
1637 auto fe_ent = DomainEleOp::getNumeredEntFiniteElementPtr()->getEnt();
1638 auto id = commonDataPtr->blocksIDmap.at(fe_ent);
1639 auto &dAta = commonDataPtr->setOfBlocksData.at(id);
1640 auto &mgis_bv = *dAta.mGisBehaviour;
1641
1642 int &size_of_vars = dAta.sizeIntVar;
1643 int &size_of_grad = dAta.sizeGradVar;
1644 int &size_of_stress = dAta.sizeStressVar;
1645
1646 auto get_tag = [&](std::string name, size_t size) {
1647 std::array<double, 9> def;
1648 std::fill(def.begin(), def.end(), 0);
1649 Tag th;
1650 CHKERR internalVarMesh.tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE,
1651 th, MB_TAG_CREAT | MB_TAG_SPARSE,
1652 def.data());
1653 return th;
1654 };
1655
1656 auto t_stress = getFTensor2FromMat<3, 3>(*(commonDataPtr->mStressPtr));
1657
1658 MatrixDouble3by3 mat(3, 3);
1659
1660 auto th_disp = get_tag(fieldName, 3);
1661 auto th_stress = get_tag(mgis_bv.thermodynamic_forces[0].name, 9);
1662 auto th_grad = get_tag(mgis_bv.gradients[0].name, 9);
1663
1664 size_t nb_gauss_pts = commonDataPtr->mGradPtr->size2();
1665 auto t_grad = getFTensor2FromMat<3, 3>(*(commonDataPtr->mGradPtr));
1666 auto t_disp = getFTensor1FromMat<3>(*(commonDataPtr->mDispPtr));
1667 CHKERR commonDataPtr->getInternalVar(fe_ent, nb_gauss_pts, size_of_vars,
1668 size_of_grad, size_of_stress);
1669
1670 MatrixDouble &mat_int = *commonDataPtr->internalVariablePtr;
1671 vector<Tag> tags_vec;
1672 bool is_large_strain = dAta.isFiniteStrain;
1673
1674 for (auto c : mgis_bv.isvs) {
1675 auto vsize = getVariableSize(c, mgis_bv.hypothesis);
1676 const size_t parav_siz = get_paraview_size(vsize);
1677 tags_vec.emplace_back(get_tag(c.name, parav_siz));
1678 }
1679
1680 if (!(side == 0 && type == MBVERTEX))
1681 MoFEMFunctionReturnHot(0);
1682
1683 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1684
1685 double coords[] = {0, 0, 0};
1686 EntityHandle vertex;
1687 for (int dd = 0; dd != 3; dd++)
1688 coords[dd] = DomainEleOp::getCoordsAtGaussPts()(gg, dd);
1689
1690 CHKERR internalVarMesh.create_vertex(coords, vertex);
1691 VectorDouble disps({t_disp(0), t_disp(1), t_disp(2)});
1692
1693 auto it = tags_vec.begin();
1694 for (auto c : mgis_bv.isvs) {
1695 auto vsize = getVariableSize(c, mgis_bv.hypothesis);
1696 const size_t parav_siz = get_paraview_size(vsize);
1697 const auto offset =
1698 getVariableOffset(mgis_bv.isvs, c.name, mgis_bv.hypothesis);
1699 auto vec =
1700 getVectorAdaptor(&mat_int.data()[gg * size_of_vars], size_of_vars);
1701 VectorDouble tag_vec = getVectorAdaptor(&vec[offset], vsize);
1702 tag_vec.resize(parav_siz);
1703
1704 CHKERR internalVarMesh.tag_set_data(*it, &vertex, 1, &*tag_vec.begin());
1705
1706 it++;
1707 }
1708
1709 // keep the convention consistent for postprocessing!
1710 array<double, 9> my_stress_vec{
1711 t_stress(0, 0), t_stress(1, 1), t_stress(2, 2),
1712 t_stress(0, 1), t_stress(1, 0), t_stress(0, 2),
1713 t_stress(2, 0), t_stress(1, 2), t_stress(2, 1)};
1714
1715 FTensor::Index<'i', 3> i;
1716 FTensor::Index<'j', 3> j;
1717
1718 array<double, 9> grad1_vec;
1719 if (is_large_strain) {
1720 FTensor::Tensor2<double, 3, 3> t_strain;
1721 t_strain(i, j) = t_grad(i, j) + kronecker_delta(i, j);
1722 grad1_vec = getVoigtVec<3>(t_strain);
1723 } else {
1724 FTensor::Tensor2_symmetric<double, 3> t_strain;
1725 t_strain(i, j) = (t_grad(i, j) || t_grad(j, i)) / 2;
1726 grad1_vec = getVoigtVecSymm<3>(t_strain);
1727 }
1728
1729 CHKERR internalVarMesh.tag_set_data(th_stress, &vertex, 1,
1730 my_stress_vec.data());
1731 CHKERR internalVarMesh.tag_set_data(th_grad, &vertex, 1,
1732 grad1_vec.data());
1733 CHKERR internalVarMesh.tag_set_data(th_disp, &vertex, 1, &*disps.begin());
1734
1735 ++t_grad;
1736 ++t_stress;
1737 ++t_disp;
1738 }
1739
1740 MoFEMFunctionReturn(0);
1741 }
1742
1743private:
1744 boost::shared_ptr<MFrontInterface::CommonData> commonDataPtr;
1745 moab::Interface &internalVarMesh;
1746 std::string fieldName;
1747};
1748
1749template <ModelHypothesis MH, AssemblyType AT>
1750MoFEMErrorCode MFrontInterfaceImpl<MH, AT>::setUpdateInternalVariablesOperators(
1751 ForcesAndSourcesCore::RuleHookFun rule, std::string field_name) {
1752 MoFEMFunctionBegin;
1753
1754 updateIntVariablesElePtr = boost::make_shared<DomainEle>(mField);
1755
1756 updateIntVariablesElePtr->getRuleHook = rule;
1757
1758 CHKERR AddHOOps<DIM, DIM, DIM>::add(
1759 updateIntVariablesElePtr->getOpPtrVector(), {H1});
1760
1761 updateIntVariablesElePtr->getOpPtrVector().push_back(
1762 new OpCalculateVectorFieldGradient<DIM, DIM>(field_name,
1763 commonDataPtr->mGradPtr));
1764 updateIntVariablesElePtr->getOpPtrVector().push_back(
1765 new OpCalculateVectorFieldValues<DIM>(field_name,
1766 commonDataPtr->mDispPtr));
1767 if (isFiniteKinematics)
1768 updateIntVariablesElePtr->getOpPtrVector().push_back(
1769 new OpUpdateVariablesFiniteStrains<MH>(field_name, commonDataPtr,
1770 monitorPtr));
1771 else
1772 updateIntVariablesElePtr->getOpPtrVector().push_back(
1773 new OpUpdateVariablesSmallStrains<MH>(field_name, commonDataPtr,
1774 monitorPtr));
1775 if (saveGauss && (MH == TRIDIMENSIONAL)) {
1776 auto &moab_gauss = *moabGaussIntPtr;
1777 updateIntVariablesElePtr->getOpPtrVector().push_back(
1778 new OpSaveGaussPts<MH>(field_name, moab_gauss, commonDataPtr));
1779 }
1780
1781 MoFEMFunctionReturn(0);
1782}
1783
1784template <bool IS_LARGE_STRAIN, ModelHypothesis MH>
1785struct OpSaveStress : public MFrontEleType<MH>::DomainEleOp {
1786 static constexpr int DIM = MFrontEleType<MH>::SPACE_DIM;
1787 using DomainEleOp = typename MFrontEleType<MH>::DomainEleOp;
1788
1789 OpSaveStress(const std::string field_name,
1790 boost::shared_ptr<MFrontInterface::CommonData> common_data_ptr)
1791 : DomainEleOp(field_name, DomainEleOp::OPROW),
1792 commonDataPtr(common_data_ptr) {
1793 std::fill(&DomainEleOp::doEntities[MBEDGE],
1794 &DomainEleOp::doEntities[MBMAXTYPE], false);
1795 }
1796 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1797 MoFEMFunctionBegin;
1798
1799 const size_t nb_gauss_pts = commonDataPtr->mGradPtr->size2();
1800 auto fe_ent = DomainEleOp::getNumeredEntFiniteElementPtr()->getEnt();
1801 auto id = commonDataPtr->blocksIDmap.at(fe_ent);
1802 auto &dAta = commonDataPtr->setOfBlocksData.at(id);
1803
1804 int &size_of_stress = dAta.sizeStressVar;
1805 int &size_of_grad = dAta.sizeGradVar;
1806
1807 CHKERR commonDataPtr->getInternalVar(fe_ent, nb_gauss_pts, dAta.sizeIntVar,
1808 dAta.sizeGradVar, dAta.sizeStressVar,
1809 IS_LARGE_STRAIN);
1810
1811 MatrixDouble &mat_stress = *commonDataPtr->mPrevStressPtr;
1812 MatrixDouble &mat_grad = *commonDataPtr->mPrevGradPtr;
1813
1814 commonDataPtr->mFullStressPtr->resize(3 * 3, nb_gauss_pts);
1815 commonDataPtr->mFullStressPtr->clear();
1816 auto t_full_stress =
1817 getFTensor2FromMat<3, 3>(*(commonDataPtr->mFullStressPtr));
1818
1819 commonDataPtr->mFullStrainPtr->resize(3 * 3, nb_gauss_pts);
1820 commonDataPtr->mFullStrainPtr->clear();
1821 auto t_full_strain =
1822 getFTensor2FromMat<3, 3>(*(commonDataPtr->mFullStrainPtr));
1823
1824 Index<'i', 3> i;
1825 Index<'j', 3> j;
1826
1827 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1828
1829 auto stress_vec = getVectorAdaptor(
1830 &mat_stress.data()[gg * size_of_stress], size_of_stress);
1831 auto grad_vec =
1832 getVectorAdaptor(&mat_grad.data()[gg * size_of_grad], size_of_grad);
1833
1834 Tensor2<double, 3, 3> t_stress;
1835 Tensor2<double, 3, 3> t_grad;
1836
1837 if constexpr (IS_LARGE_STRAIN) {
1838 if constexpr (DIM == 3) {
1839 t_stress = Tensor2<double, 3, 3>(VOIGT_VEC_3D(stress_vec));
1840 t_grad = Tensor2<double, 3, 3>(VOIGT_VEC_3D(grad_vec));
1841 } else if constexpr (DIM == 2) {
1842 t_stress = Tensor2<double, 3, 3>(VOIGT_VEC_2D_FULL(stress_vec));
1843 t_grad = Tensor2<double, 3, 3>(VOIGT_VEC_2D_FULL(grad_vec));
1844 }
1845 } else {
1846 if constexpr (DIM == 3) {
1847 t_stress = to_non_symm(
1848 Tensor2_symmetric<double, 3>(VOIGT_VEC_SYMM_3D(stress_vec)));
1849 t_grad = to_non_symm(
1850 Tensor2_symmetric<double, 3>(VOIGT_VEC_SYMM_3D(grad_vec)));
1851 } else if constexpr (DIM == 2) {
1852 t_stress = to_non_symm(
1853 Tensor2_symmetric<double, 3>(VOIGT_VEC_SYMM_2D_FULL(stress_vec)));
1854 t_grad = to_non_symm(
1855 Tensor2_symmetric<double, 3>(VOIGT_VEC_SYMM_2D_FULL(grad_vec)));
1856 }
1857 }
1858
1859 t_full_stress(i, j) = t_stress(i, j);
1860 t_full_strain(i, j) = t_grad(i, j);
1861
1862 ++t_full_stress;
1863 ++t_full_strain;
1864 }
1865
1866 MoFEMFunctionReturn(0);
1867 }
1868
1869private:
1870 boost::shared_ptr<MFrontInterface::CommonData> commonDataPtr;
1871};
1872
1873template <ModelHypothesis MH, AssemblyType AT>
1874MoFEMErrorCode MFrontInterfaceImpl<MH, AT>::setPostProcessOperators(
1875 ForcesAndSourcesCore::RuleHookFun rule, std::string fe_name,
1876 std::string field_name, int order) {
1877 MoFEMFunctionBegin;
1878
1879 postProcFe = boost::make_shared<PostProcDomainEle>(mField);
1880
1881 auto &pip = postProcFe->getOpPtrVector();
1882
1883 pip.push_back(new OpCalculateVectorFieldValues<DIM>(field_name,
1884 commonDataPtr->mDispPtr));
1885
1886 auto entity_data_l2 = boost::make_shared<EntitiesFieldData>(MBENTITYSET);
1887 auto mass_ptr = boost::make_shared<MatrixDouble>();
1888 auto strain_coeffs_ptr = boost::make_shared<MatrixDouble>();
1889 auto stress_coeffs_ptr = boost::make_shared<MatrixDouble>();
1890
1891 auto op_this = new OpLoopThis<DomainEle>(mField, fe_name, Sev::noisy);
1892 pip.push_back(op_this);
1893 pip.push_back(new OpDGProjectionEvaluation(commonDataPtr->mFullStrainPtr,
1894 strain_coeffs_ptr, entity_data_l2,
1895 fieldBase, L2));
1896 pip.push_back(new OpDGProjectionEvaluation(commonDataPtr->mFullStressPtr,
1897 stress_coeffs_ptr, entity_data_l2,
1898 fieldBase, L2));
1899
1900 auto fe_physics_ptr = op_this->getThisFEPtr();
1901 fe_physics_ptr->getRuleHook = rule;
1902
1903 fe_physics_ptr->getOpPtrVector().push_back(new OpDGProjectionMassMatrix(
1904 order, mass_ptr, entity_data_l2, fieldBase, L2));
1905 if (isFiniteKinematics) {
1906 fe_physics_ptr->getOpPtrVector().push_back(
1907 new OpSaveStress<true, MH>(field_name, commonDataPtr));
1908 } else {
1909 fe_physics_ptr->getOpPtrVector().push_back(
1910 new OpSaveStress<false, MH>(field_name, commonDataPtr));
1911 }
1912 fe_physics_ptr->getOpPtrVector().push_back(new OpDGProjectionCoefficients(
1913 commonDataPtr->mFullStrainPtr, strain_coeffs_ptr, mass_ptr,
1914 entity_data_l2, fieldBase, L2));
1915 fe_physics_ptr->getOpPtrVector().push_back(new OpDGProjectionCoefficients(
1916 commonDataPtr->mFullStressPtr, stress_coeffs_ptr, mass_ptr,
1917 entity_data_l2, fieldBase, L2));
1918
1919 using OpPPMapVec = OpPostProcMapInMoab<DIM, DIM>;
1920 using OpPPMapTen = OpPostProcMapInMoab<3, 3>;
1921
1922 pip.push_back(new OpPPMapVec(
1923
1924 postProcFe->getPostProcMesh(), postProcFe->getMapGaussPts(),
1925
1926 {},
1927
1928 {{field_name, commonDataPtr->mDispPtr}},
1929
1930 {},
1931
1932 {}
1933
1934 ));
1935
1936 pip.push_back(new OpPPMapTen(
1937
1938 postProcFe->getPostProcMesh(), postProcFe->getMapGaussPts(),
1939
1940 {},
1941
1942 {},
1943
1944 {{"STRAIN", commonDataPtr->mFullStrainPtr},
1945 {"STRESS", commonDataPtr->mFullStressPtr}},
1946
1947 {}
1948
1949 ));
1950
1951 MoFEMFunctionReturn(0);
1952}
1953
1954template <ModelHypothesis MH, AssemblyType AT>
1955MoFEMErrorCode MFrontInterfaceImpl<MH, AT>::postProcess(int step,
1956 SmartPetscObj<DM> dm,
1957 string fe_name) {
1958 MoFEMFunctionBegin;
1959
1960 auto make_vtks = [&]() {
1961 MoFEMFunctionBegin;
1962
1963 if (saveDomain) {
1964 CHKERR DMoFEMLoopFiniteElements(dm, fe_name, postProcFe);
1965 CHKERR postProcFe->writeFile("out_" + optionsPrefix +
1966 boost::lexical_cast<std::string>(step) +
1967 ".h5m");
1968 }
1969
1970 if (saveGauss) {
1971 string file_name = "out_" + optionsPrefix + "gauss_" +
1972 boost::lexical_cast<std::string>(step) + ".h5m";
1973
1974 CHKERR moabGaussIntPtr->write_file(file_name.c_str(), "MOAB",
1975 "PARALLEL=WRITE_PART");
1976 CHKERR moabGaussIntPtr->delete_mesh();
1977 }
1978
1979 MoFEMFunctionReturn(0);
1980 };
1981
1982 CHKERR make_vtks();
1983
1984 MoFEMFunctionReturn(0);
1985};
1986
1987template <ModelHypothesis MH, AssemblyType AT>
1988MoFEMErrorCode
1989MFrontInterfaceImpl<MH, AT>::updateInternalVariables(SmartPetscObj<DM> dm,
1990 std::string fe_name) {
1991 MoFEMFunctionBegin;
1992 CHKERR DMoFEMLoopFiniteElements(dm, fe_name, updateIntVariablesElePtr);
1993 MoFEMFunctionReturn(0);
1994}
1995
1996} // namespace MoFEM
1997
1998#else // WITH_MGIS
1999
2000namespace MoFEM {
2001
2005
2006template <ModelHypothesis MH, AssemblyType AT>
2007MFrontInterfaceImpl<MH, AT>::MFrontInterfaceImpl(MoFEM::Interface &m_field)
2008 : mField(m_field) {}
2009
2010template <ModelHypothesis MH, AssemblyType AT>
2011MoFEMErrorCode MFrontInterfaceImpl<MH, AT>::getCommandLineParameters() {
2012 MoFEMFunctionBegin;
2013 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_INSTALLED,
2014 "MFront has not been enabled in this build of MoFEM");
2015 MoFEMFunctionReturn(0);
2016}
2017
2018template <ModelHypothesis MH, AssemblyType AT>
2019MoFEMErrorCode MFrontInterfaceImpl<MH, AT>::opFactoryDomainRhs(
2020 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
2021 std::string field_name) {
2022 MoFEMFunctionBegin;
2023 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_INSTALLED,
2024 "MFront has not been enabled in this build of MoFEM");
2025 MoFEMFunctionReturn(0);
2026}
2027
2028template <ModelHypothesis MH, AssemblyType AT>
2029MoFEMErrorCode MFrontInterfaceImpl<MH, AT>::opFactoryDomainLhs(
2030 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
2031 std::string field_name) {
2032 MoFEMFunctionBegin;
2033 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_INSTALLED,
2034 "MFront has not been enabled in this build of MoFEM");
2035 MoFEMFunctionReturn(0);
2036}
2037
2038template <ModelHypothesis MH, AssemblyType AT>
2039MoFEMErrorCode MFrontInterfaceImpl<MH, AT>::setUpdateInternalVariablesOperators(
2040 ForcesAndSourcesCore::RuleHookFun rule, std::string field_name) {
2041 MoFEMFunctionBegin;
2042 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_INSTALLED,
2043 "MFront has not been enabled in this build of MoFEM");
2044 MoFEMFunctionReturn(0);
2045}
2046
2047template <ModelHypothesis MH, AssemblyType AT>
2048MoFEMErrorCode MFrontInterfaceImpl<MH, AT>::setPostProcessOperators(
2049 ForcesAndSourcesCore::RuleHookFun rule, std::string fe_name,
2050 std::string field_name, int order) {
2051 MoFEMFunctionBegin;
2052 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_INSTALLED,
2053 "MFront has not been enabled in this build of MoFEM");
2054 MoFEMFunctionReturn(0);
2055}
2056
2057template <ModelHypothesis MH, AssemblyType AT>
2058MoFEMErrorCode MFrontInterfaceImpl<MH, AT>::postProcess(int step,
2059 SmartPetscObj<DM> dm,
2060 string fe_name) {
2061 MoFEMFunctionBegin;
2062 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_INSTALLED,
2063 "MFront has not been enabled in this build of MoFEM");
2064 MoFEMFunctionReturn(0);
2065};
2066
2067template <ModelHypothesis MH, AssemblyType AT>
2068MoFEMErrorCode
2069MFrontInterfaceImpl<MH, AT>::updateInternalVariables(SmartPetscObj<DM> dm,
2070 std::string fe_name) {
2071 MoFEMFunctionBegin;
2072 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_INSTALLED,
2073 "MFront has not been enabled in this build of MoFEM");
2074 MoFEMFunctionReturn(0);
2075}
2076
2077} // namespace MoFEM
2078
2079#endif // WITH_MGIS
2080
2081namespace MoFEM {
2082
2083template struct MFrontInterfaceImpl<TRIDIMENSIONAL, AssemblyType::PETSC>;
2084template struct MFrontInterfaceImpl<AXISYMMETRICAL, AssemblyType::PETSC>;
2085template struct MFrontInterfaceImpl<PLANESTRAIN, AssemblyType::PETSC>;
2086
2087template struct MFrontInterfaceImpl<TRIDIMENSIONAL, AssemblyType::BLOCK_SCHUR>;
2088template struct MFrontInterfaceImpl<AXISYMMETRICAL, AssemblyType::BLOCK_SCHUR>;
2089template struct MFrontInterfaceImpl<PLANESTRAIN, AssemblyType::BLOCK_SCHUR>;
2090
2091} // namespace MoFEM
setBlocks
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData > > &block_sets_ptr)
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
is_large_strain
PetscBool is_large_strain
Definition contact.cpp:92
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_NOT_INSTALLED
@ MOFEM_NOT_INSTALLED
Definition definitions.h:37
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:19
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:39
D
double D
Definition initial_diffusion.cpp:44
n
const double n
refractive index of diffusive medium
Definition initial_diffusion.cpp:38
J
FTensor::Index< 'J', DIM1 > J
Definition level_set.cpp:30
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:19
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::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition Types.hpp:77
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition Types.hpp:68
MoFEM::Types::MatrixAdaptor
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Matrix adaptor.
Definition Types.hpp:132
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:31
MoFEM::type_from_handle
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition Templates.hpp:1992
MoFEM::getVoigtVecAxisymm
auto getVoigtVecAxisymm(T &t_mat, const double hoop_term)
Definition Templates.hpp:1425
MoFEM::to_non_symm
auto to_non_symm(const FTensor::Tensor2_symmetric< T, DIM > &symm)
Definition Templates.hpp:1444
MoFEM::createMFrontInterface
boost::shared_ptr< MFrontInterface > createMFrontInterface(MoFEM::Interface &m_field, ModelHypothesis mh, AssemblyType at=AssemblyType::PETSC)
create mfront interface
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:1433
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:500
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:547
sdf.r
int r
Definition sdf.py:8
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:60
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:152
MoFEM::EntitiesFieldData::EntData::getFTensor1DiffN
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition EntitiesFieldData.cpp:526
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition EntitiesFieldData.hpp:1692
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:1486
MoFEM::ForcesAndSourcesCore::RuleHookFun
boost::function< int(int order_row, int order_col, int order_data)> RuleHookFun
Definition ForcesAndSourcesCore.hpp:79
MoFEM::MFrontInterface::CommonData::CommonData
CommonData(MoFEM::Interface &m_field)
Definition MFrontInterface.cpp:2003
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 &row_data, EntData &col_data)
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
Generated by Doxygen 1.9.8 and hosted at University of Glasgow