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