v0.12.1
plastic.cpp
Go to the documentation of this file.
1 /**
2  * \file plastic.cpp
3  * \example plastic.cpp
4  *
5  * Plasticity in 2d and 3d
6  *
7  */
8 
9 /* This file is part of MoFEM.
10  * MoFEM is free software: you can redistribute it and/or modify it under
11  * the terms of the GNU Lesser General Public License as published by the
12  * Free Software Foundation, either version 3 of the License, or (at your
13  * option) any later version.
14  *
15  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
16  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18  * License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public
21  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
22 
23 #ifndef EXECUTABLE_DIMENSION
24 #define EXECUTABLE_DIMENSION 3
25 #endif
26 
27 #include <MoFEM.hpp>
28 #include <MatrixFunction.hpp>
29 #include <IntegrationRules.hpp>
30 
31 using namespace MoFEM;
32 
33 template <int DIM> struct ElementsAndOps {};
34 
35 template <> struct ElementsAndOps<2> {
41 };
42 
43 template <> struct ElementsAndOps<3> {
49 };
50 
51 constexpr int SPACE_DIM =
52  EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
53 
54 constexpr EntityType boundary_ent = SPACE_DIM == 3 ? MBTRI : MBEDGE;
61 
64 
65 //! [Body force]
67  GAUSS>::OpSource<1, SPACE_DIM>;
68 //! [Body force]
69 
70 //! [Only used with Hooke equation (linear material model)]
72  GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
75 //! [Only used with Hooke equation (linear material model)]
76 
77 //! [Only used for dynamics]
82 //! [Only used for dynamics]
83 
84 //! [Only used with Hencky/nonlinear material]
86  GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
89 //! [Only used with Hencky/nonlinear material]
90 
91 //! [Essential boundary conditions]
98 //! [Essential boundary conditions]
100 
101 PetscBool is_large_strains = PETSC_TRUE;
102 
103 double scale = 1.;
104 
105 double young_modulus = 206913;
106 double poisson_ratio = 0.29;
107 double rho = 0;
108 double sigmaY = 450;
109 double H = 129;
110 double visH = 0;
111 double cn = 1;
112 double Qinf = 265;
113 double b_iso = 16.93;
114 int order = 2;
115 
116 inline long double hardening(long double tau, double temp) {
117  return H * tau + Qinf * (1. - std::exp(-b_iso * tau)) + sigmaY;
118 }
119 
120 inline long double hardening_dtau(long double tau, double temp) {
121  return H + Qinf * b_iso * std::exp(-b_iso * tau);
122 }
123 
124 #include <HenckyOps.hpp>
125 #include <PlasticOps.hpp>
126 #include <OpPostProcElastic.hpp>
127 
128 using namespace PlasticOps;
129 using namespace HenckyOps;
130 struct Example {
131 
132  Example(MoFEM::Interface &m_field) : mField(m_field) {}
133 
135 
136 private:
137  MoFEM::Interface &mField;
138 
141  MoFEMErrorCode bC();
142  MoFEMErrorCode OPs();
143  MoFEMErrorCode tsSolve();
144 
145  boost::shared_ptr<PlasticOps::CommonData> commonPlasticDataPtr;
146  boost::shared_ptr<HenckyOps::CommonData> commonHenckyDataPtr;
147  boost::shared_ptr<PostProcEle> postProcFe;
148  boost::shared_ptr<DomainEle> reactionFe;
149  boost::shared_ptr<DomainEle> feAxiatorRhs;
150  boost::shared_ptr<DomainEle> feAxiatorLhs;
151 
152  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
153  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
154  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
155 
156  boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
157  boost::shared_ptr<std::vector<unsigned char>> reactionMarker;
158 
159  std::vector<FTensor::Tensor1<double, 3>> bodyForces;
160 };
161 
162 //! [Run problem]
165  CHKERR setupProblem();
166  CHKERR createCommonData();
167  CHKERR bC();
168  CHKERR OPs();
169  CHKERR tsSolve();
171 }
172 //! [Run problem]
173 
174 //! [Set up problem]
177  Simple *simple = mField.getInterface<Simple>();
178  // Add field
180  constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
181  CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
182  CHKERR simple->addDomainField("TAU", L2, base, 1);
183  CHKERR simple->addDomainField("EP", L2, base, size_symm);
184  CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
185  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
186  CHKERR simple->setFieldOrder("U", order);
187  CHKERR simple->setFieldOrder("TAU", order - 1);
188  CHKERR simple->setFieldOrder("EP", order - 1);
189  CHKERR simple->setUp();
191 }
192 //! [Set up problem]
193 
194 //! [Create common data]
197 
198  auto get_command_line_parameters = [&]() {
200  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-scale", &scale, PETSC_NULL);
201  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-rho", &rho, PETSC_NULL);
202  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
203  &young_modulus, PETSC_NULL);
204  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
205  &poisson_ratio, PETSC_NULL);
206  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening", &H, PETSC_NULL);
207  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening_viscous", &visH,
208  PETSC_NULL);
209  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-yield_stress", &sigmaY,
210  PETSC_NULL);
211  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn", &cn, PETSC_NULL);
212  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-Qinf", &Qinf, PETSC_NULL);
213  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-b_iso", &b_iso, PETSC_NULL);
214 
215  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-large_strains",
216  &is_large_strains, PETSC_NULL);
217 
218  MOFEM_LOG("EXAMPLE", Sev::inform) << "Young modulus " << young_modulus;
219  MOFEM_LOG("EXAMPLE", Sev::inform) << "Poisson ratio " << poisson_ratio;
220  MOFEM_LOG("EXAMPLE", Sev::inform) << "Yield stress " << sigmaY;
221  MOFEM_LOG("EXAMPLE", Sev::inform) << "Hardening " << H;
222  MOFEM_LOG("EXAMPLE", Sev::inform) << "Viscous hardening " << visH;
223  MOFEM_LOG("EXAMPLE", Sev::inform) << "Saturation yield stress " << Qinf;
224  MOFEM_LOG("EXAMPLE", Sev::inform) << "Saturation exponent " << b_iso;
225  MOFEM_LOG("EXAMPLE", Sev::inform) << "cn " << cn;
226 
227  PetscBool is_scale = PETSC_TRUE;
228  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_scale", &is_scale,
229  PETSC_NULL);
230  if (is_scale) {
232  young_modulus *= scale;
233  rho *= scale;
234  sigmaY *= scale;
235  H *= scale;
236  Qinf *= scale;
237  visH *= scale;
238 
239  MOFEM_LOG("EXAMPLE", Sev::inform)
240  << "Scaled Young modulus " << young_modulus;
241  MOFEM_LOG("EXAMPLE", Sev::inform)
242  << "Scaled Poisson ratio " << poisson_ratio;
243  MOFEM_LOG("EXAMPLE", Sev::inform) << "Scaled Yield stress " << sigmaY;
244  MOFEM_LOG("EXAMPLE", Sev::inform) << "Scaled Hardening " << H;
245  MOFEM_LOG("EXAMPLE", Sev::inform) << "Scaled Viscous hardening " << visH;
246  MOFEM_LOG("EXAMPLE", Sev::inform)
247  << "Scaled Saturation yield stress " << Qinf;
248  }
249 
251  };
252 
253  auto set_matrial_stiffness = [&]() {
260  const double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
261  const double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
262 
263  // Plane stress or when 1, plane strain or 3d
264  const double A = (SPACE_DIM == 2)
265  ? 2 * shear_modulus_G /
266  (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
267  : 1;
268 
269  auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
270  *commonPlasticDataPtr->mDPtr);
271  auto t_D_axiator = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
272  *commonPlasticDataPtr->mDPtr_Axiator);
273  auto t_D_deviator = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
274  *commonPlasticDataPtr->mDPtr_Deviator);
275 
276  constexpr double third = boost::math::constants::third<double>();
277  t_D_axiator(i, j, k, l) = A *
278  (bulk_modulus_K - (2. / 3.) * shear_modulus_G) *
279  t_kd(i, j) * t_kd(k, l);
280  t_D_deviator(i, j, k, l) =
281  2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
282  t_D(i, j, k, l) = t_D_axiator(i, j, k, l) + t_D_deviator(i, j, k, l);
283 
285  };
286 
287  commonPlasticDataPtr = boost::make_shared<PlasticOps::CommonData>();
288  constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
289  commonPlasticDataPtr->mDPtr = boost::make_shared<MatrixDouble>();
290  commonPlasticDataPtr->mDPtr->resize(size_symm * size_symm, 1);
291  commonPlasticDataPtr->mDPtr_Axiator = boost::make_shared<MatrixDouble>();
292  commonPlasticDataPtr->mDPtr_Axiator->resize(size_symm * size_symm, 1);
293  commonPlasticDataPtr->mDPtr_Deviator = boost::make_shared<MatrixDouble>();
294  commonPlasticDataPtr->mDPtr_Deviator->resize(size_symm * size_symm, 1);
295 
296  commonPlasticDataPtr->mGradPtr = boost::make_shared<MatrixDouble>();
297  commonPlasticDataPtr->mStrainPtr = boost::make_shared<MatrixDouble>();
298  commonPlasticDataPtr->mStressPtr = boost::make_shared<MatrixDouble>();
299 
300  CHKERR get_command_line_parameters();
301  CHKERR set_matrial_stiffness();
302 
303  if (is_large_strains) {
304  commonHenckyDataPtr = boost::make_shared<HenckyOps::CommonData>();
305  commonHenckyDataPtr->matGradPtr = commonPlasticDataPtr->mGradPtr;
306  commonHenckyDataPtr->matDPtr = commonPlasticDataPtr->mDPtr;
307  commonHenckyDataPtr->matLogCPlastic =
308  commonPlasticDataPtr->getPlasticStrainPtr();
309  commonPlasticDataPtr->mStrainPtr = commonHenckyDataPtr->getMatLogC();
310  commonPlasticDataPtr->mStressPtr =
311  commonHenckyDataPtr->getMatHenckyStress();
312  }
313 
315 }
316 //! [Create common data]
317 
318 //! [Boundary condition]
321 
322  auto simple = mField.getInterface<Simple>();
323  auto bc_mng = mField.getInterface<BcManager>();
324  auto prb_mng = mField.getInterface<ProblemsManager>();
325 
326  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
327  "U", 0, 0);
328  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
329  "U", 1, 1);
330  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
331  "U", 2, 2);
332  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
333  "REMOVE_ALL", "U", 0, 3);
334 
335  CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_X", "U",
336  0, 0);
337  CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_Y", "U",
338  1, 1);
339  CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_Z", "U",
340  2, 2);
341  CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
342  "U", 0, 3);
343 
344  auto &bc_map = bc_mng->getBcMapByBlockName();
345  boundaryMarker = bc_mng->getMergedBlocksMarker(vector<string>{"FIX_"});
346 
347  CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "REACTION",
348  "U", 0, 3);
349 
350  for (auto bc : bc_map)
351  MOFEM_LOG("EXAMPLE", Sev::verbose) << "Marker " << bc.first;
352 
353  // OK. We have problem with GMesh, it adding empty characters at the end of
354  // block. So first block is search by regexp. popMarkDOFsOnEntities should
355  // work with regexp.
356  std::string reaction_block_set;
357  for (auto bc : bc_map) {
358  if (bc_mng->checkBlock(bc, "REACTION")) {
359  reaction_block_set = bc.first;
360  break;
361  }
362  }
363 
364  if (auto bc = bc_mng->popMarkDOFsOnEntities(reaction_block_set)) {
365  reactionMarker = bc->getBcMarkersPtr();
366 
367  // Only take reaction from nodes
368  Range nodes;
369  CHKERR mField.get_moab().get_entities_by_type(0, MBVERTEX, nodes, true);
370  CHKERR prb_mng->markDofs(simple->getProblemName(), ROW,
371  ProblemsManager::MarkOP::AND, nodes,
372  *reactionMarker);
373 
374  } else {
375  MOFEM_LOG("EXAMPLE", Sev::warning) << "REACTION blockset does not exist";
376  }
377 
378  if (!reactionMarker) {
379  MOFEM_LOG("EXAMPLE", Sev::warning) << "REACTION blockset does not exist";
380  }
381 
383 }
384 //! [Boundary condition]
385 
386 //! [Push operators to pipeline]
389  auto pipeline_mng = mField.getInterface<PipelineManager>();
390  auto simple = mField.getInterface<Simple>();
391  auto bc_mng = mField.getInterface<BcManager>();
392 
393  feAxiatorLhs = boost::make_shared<DomainEle>(mField);
394  feAxiatorRhs = boost::make_shared<DomainEle>(mField);
395  auto integration_rule_axiator = [](int, int, int approx_order) {
396  return 2 * (approx_order - 1);
397  };
398  feAxiatorLhs->getRuleHook = integration_rule_axiator;
399  feAxiatorRhs->getRuleHook = integration_rule_axiator;
400 
401  auto integration_rule_deviator = [](int o_row, int o_col, int approx_order) {
402  return 2 * (approx_order - 1);
403  };
404  auto integration_rule_bc = [](int, int, int approx_order) {
405  return 2 * approx_order;
406  };
407 
408  auto add_domain_base_ops = [&](auto &pipeline) {
410 
411  if (SPACE_DIM == 2) {
412  auto det_ptr = boost::make_shared<VectorDouble>();
413  auto jac_ptr = boost::make_shared<MatrixDouble>();
414  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
415  pipeline.push_back(new OpCalculateHOJacForFace(jac_ptr));
416  pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
417  pipeline.push_back(new OpSetInvJacH1ForFace(inv_jac_ptr));
418  }
419 
420  pipeline.push_back(new OpCalculateScalarFieldValuesDot(
421  "TAU", commonPlasticDataPtr->getPlasticTauDotPtr()));
423  "EP", commonPlasticDataPtr->getPlasticStrainPtr()));
425  "EP", commonPlasticDataPtr->getPlasticStrainDotPtr()));
427  "U", commonPlasticDataPtr->mGradPtr));
428  pipeline.push_back(new OpCalculateScalarFieldValues(
429  "TAU", commonPlasticDataPtr->getPlasticTauPtr()));
430 
432  };
433 
434  auto add_domain_stress_ops = [&](auto &pipeline, auto m_D_ptr) {
436 
437  if (is_large_strains) {
438 
439  if (commonPlasticDataPtr->mGradPtr != commonHenckyDataPtr->matGradPtr)
440  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
441  "Wrong pointer for grad");
442 
443  pipeline.push_back(
444  new OpCalculateEigenVals<SPACE_DIM>("U", commonHenckyDataPtr));
445  pipeline.push_back(
446  new OpCalculateLogC<SPACE_DIM>("U", commonHenckyDataPtr));
447  pipeline.push_back(
448  new OpCalculateLogC_dC<SPACE_DIM>("U", commonHenckyDataPtr));
449  pipeline.push_back(new OpCalculateHenckyPlasticStress<SPACE_DIM>(
450  "U", commonHenckyDataPtr, m_D_ptr));
451  pipeline.push_back(
452  new OpCalculatePiolaStress<SPACE_DIM>("U", commonHenckyDataPtr));
453 
454  } else {
455  pipeline.push_back(
456  new OpSymmetrizeTensor<SPACE_DIM>("U", commonPlasticDataPtr->mGradPtr,
457  commonPlasticDataPtr->mStrainPtr));
458  pipeline.push_back(
459  new OpPlasticStress("U", commonPlasticDataPtr, m_D_ptr, 1));
460  }
461 
462  if (m_D_ptr != commonPlasticDataPtr->mDPtr_Axiator)
463  pipeline.push_back(
464  new OpCalculatePlasticSurface("U", commonPlasticDataPtr));
465 
467  };
468 
469  auto add_domain_ops_lhs_mechanical = [&](auto &pipeline, auto m_D_ptr) {
471  pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
472 
473  if (is_large_strains) {
474  pipeline.push_back(
475  new OpHenckyTangent<SPACE_DIM>("U", commonHenckyDataPtr, m_D_ptr));
476  pipeline.push_back(
477  new OpKPiola("U", "U", commonHenckyDataPtr->getMatTangent()));
479  "U", "EP", commonPlasticDataPtr, commonHenckyDataPtr, m_D_ptr));
480  } else {
481  pipeline.push_back(new OpKCauchy("U", "U", m_D_ptr));
482  pipeline.push_back(new OpCalculatePlasticInternalForceLhs_dEP(
483  "U", "EP", commonPlasticDataPtr, m_D_ptr));
484  }
485 
486  pipeline.push_back(new OpUnSetBc("U"));
488  };
489 
490  auto add_domain_ops_rhs_mechanical = [&](auto &pipeline) {
492  pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
493 
495  const std::string block_name = "BODY_FORCE";
496  if (it->getName().compare(0, block_name.size(), block_name) == 0) {
497  std::vector<double> attr;
498  CHKERR it->getAttributes(attr);
499  if (attr.size() == 3) {
500  bodyForces.push_back(
501  FTensor::Tensor1<double, 3>{attr[0], attr[1], attr[2]});
502  } else {
503  SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
504  "Should be three atributes in BODYFORCE blockset, but is %d",
505  attr.size());
506  }
507  }
508  }
509 
510  auto get_body_force = [this](const double, const double, const double) {
511  auto *pipeline_mng = mField.getInterface<PipelineManager>();
514  t_source(i) = 0;
515  auto fe_domain_rhs = pipeline_mng->getDomainRhsFE();
516  const auto time = fe_domain_rhs->ts_t;
517  // hardcoded gravity load
518  for (auto &t_b : bodyForces)
519  t_source(i) += (scale * t_b(i)) * time;
520  return t_source;
521  };
522 
523  pipeline.push_back(new OpBodyForce("U", get_body_force));
524 
525  // Calculate internal forece
526  if (is_large_strains) {
527  pipeline.push_back(new OpInternalForcePiola(
528  "U", commonHenckyDataPtr->getMatFirstPiolaStress()));
529  } else {
530  pipeline.push_back(
531  new OpInternalForceCauchy("U", commonPlasticDataPtr->mStressPtr));
532  }
533 
534  pipeline.push_back(new OpUnSetBc("U"));
536  };
537 
538  auto add_domain_ops_lhs_constrain = [&](auto &pipeline, auto m_D_ptr) {
540  pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
541 
542  if (is_large_strains) {
543  pipeline.push_back(
544  new OpHenckyTangent<SPACE_DIM>("U", commonHenckyDataPtr));
545  pipeline.push_back(new OpCalculateContrainsLhs_LogStrain_dU(
546  "TAU", "U", commonPlasticDataPtr, commonHenckyDataPtr, m_D_ptr));
547  pipeline.push_back(new OpCalculatePlasticFlowLhs_LogStrain_dU(
548  "EP", "U", commonPlasticDataPtr, commonHenckyDataPtr, m_D_ptr));
549  } else {
550  pipeline.push_back(new OpCalculatePlasticFlowLhs_dU(
551  "EP", "U", commonPlasticDataPtr, m_D_ptr));
552  pipeline.push_back(new OpCalculateContrainsLhs_dU(
553  "TAU", "U", commonPlasticDataPtr, m_D_ptr));
554  }
555 
556  pipeline.push_back(new OpCalculatePlasticFlowLhs_dEP(
557  "EP", "EP", commonPlasticDataPtr, m_D_ptr));
558  pipeline.push_back(
559  new OpCalculatePlasticFlowLhs_dTAU("EP", "TAU", commonPlasticDataPtr));
560  pipeline.push_back(new OpCalculateContrainsLhs_dEP(
561  "TAU", "EP", commonPlasticDataPtr, m_D_ptr));
562  pipeline.push_back(
563  new OpCalculateContrainsLhs_dTAU("TAU", "TAU", commonPlasticDataPtr));
564 
565  pipeline.push_back(new OpUnSetBc("U"));
567  };
568 
569  auto add_domain_ops_rhs_constrain = [&](auto &pipeline) {
571  pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
572 
573 
574  pipeline.push_back(
575  new OpCalculatePlasticFlowRhs("EP", commonPlasticDataPtr));
576  pipeline.push_back(
577  new OpCalculateContrainsRhs("TAU", commonPlasticDataPtr));
578 
580  };
581 
582  auto add_boundary_ops_lhs_mechanical = [&](auto &pipeline) {
584  auto &bc_map = mField.getInterface<BcManager>()->getBcMapByBlockName();
585  for (auto bc : bc_map) {
586  if (bc_mng->checkBlock(bc, "FIX_")){
587  pipeline.push_back(
588  new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
589  pipeline.push_back(new OpBoundaryMass(
590  "U", "U", [](double, double, double) { return 1.; },
591  bc.second->getBcEntsPtr()));
592  pipeline.push_back(new OpUnSetBc("U"));
593  }
594  }
596  };
597 
598  auto add_boundary_ops_rhs_mechanical = [&](auto &pipeline) {
600 
601  auto get_time = [&](double, double, double) {
602  auto *pipeline_mng = mField.getInterface<PipelineManager>();
603  auto &fe_domain_rhs = pipeline_mng->getBoundaryRhsFE();
604  return fe_domain_rhs->ts_t;
605  };
606 
607  auto get_time_scaled = [&](double, double, double) {
608  auto *pipeline_mng = mField.getInterface<PipelineManager>();
609  auto &fe_domain_rhs = pipeline_mng->getBoundaryRhsFE();
610  return fe_domain_rhs->ts_t * scale;
611  };
612 
613  auto get_minus_time = [&](double, double, double) {
614  auto *pipeline_mng = mField.getInterface<PipelineManager>();
615  auto &fe_domain_rhs = pipeline_mng->getBoundaryRhsFE();
616  return -fe_domain_rhs->ts_t;
617  };
618 
619  auto time_scaled = [&](double, double, double) {
620  auto *pipeline_mng = mField.getInterface<PipelineManager>();
621  auto &fe_domain_rhs = pipeline_mng->getBoundaryRhsFE();
622  return -fe_domain_rhs->ts_t;
623  };
624 
625  pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
626 
628  if (it->getName().compare(0, 5, "FORCE") == 0) {
629  Range force_edges;
630  std::vector<double> attr_vec;
631  CHKERR it->getMeshsetIdEntitiesByDimension(
632  mField.get_moab(), SPACE_DIM - 1, force_edges, true);
633  it->getAttributes(attr_vec);
634  if (attr_vec.size() < SPACE_DIM)
635  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
636  "Wrong size of boundary attributes vector. Set right block "
637  "size attributes.");
638  auto force_vec_ptr = boost::make_shared<MatrixDouble>(SPACE_DIM, 1);
639  std::copy(&attr_vec[0], &attr_vec[SPACE_DIM],
640  force_vec_ptr->data().begin());
641  pipeline.push_back(
642  new OpBoundaryVec("U", force_vec_ptr, get_time_scaled,
643  boost::make_shared<Range>(force_edges)));
644  }
645  }
646 
647  pipeline.push_back(new OpUnSetBc("U"));
648 
649  auto u_mat_ptr = boost::make_shared<MatrixDouble>();
650  pipeline.push_back(
651  new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_mat_ptr));
652 
653  for (auto &bc : mField.getInterface<BcManager>()->getBcMapByBlockName()) {
654  if (bc_mng->checkBlock(bc, "FIX_")) {
655  pipeline.push_back(
656  new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
657  auto attr_vec = boost::make_shared<MatrixDouble>(SPACE_DIM, 1);
658  attr_vec->clear();
659  if (bc.second->bcAttributes.size() < SPACE_DIM)
660  SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
661  "Wrong size of boundary attributes vector. Set right block "
662  "size attributes. Size of attributes %d",
663  bc.second->bcAttributes.size());
664  std::copy(&bc.second->bcAttributes[0],
665  &bc.second->bcAttributes[SPACE_DIM],
666  attr_vec->data().begin());
667 
668  pipeline.push_back(new OpBoundaryVec("U", attr_vec, time_scaled,
669  bc.second->getBcEntsPtr()));
670  pipeline.push_back(new OpBoundaryInternal(
671  "U", u_mat_ptr, [](double, double, double) { return 1.; },
672  bc.second->getBcEntsPtr()));
673 
674  pipeline.push_back(new OpUnSetBc("U"));
675  }
676  }
677 
679  };
680 
681  // Axiator
682  CHKERR add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
683  CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainLhsPipeline(),
684  commonPlasticDataPtr->mDPtr_Deviator);
685  CHKERR add_domain_ops_lhs_mechanical(pipeline_mng->getOpDomainLhsPipeline(),
686  commonPlasticDataPtr->mDPtr_Deviator);
687  CHKERR add_domain_ops_lhs_constrain(pipeline_mng->getOpDomainLhsPipeline(),
688  commonPlasticDataPtr->mDPtr_Deviator);
689  CHKERR add_boundary_ops_lhs_mechanical(
690  pipeline_mng->getOpBoundaryLhsPipeline());
691 
692  CHKERR add_domain_base_ops(feAxiatorLhs->getOpPtrVector());
693  CHKERR add_domain_stress_ops(feAxiatorLhs->getOpPtrVector(),
694  commonPlasticDataPtr->mDPtr_Axiator);
695  CHKERR add_domain_ops_lhs_mechanical(feAxiatorLhs->getOpPtrVector(),
696  commonPlasticDataPtr->mDPtr_Axiator);
697 
698  CHKERR add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
699  CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainRhsPipeline(),
700  commonPlasticDataPtr->mDPtr_Deviator);
701  CHKERR add_domain_ops_rhs_mechanical(pipeline_mng->getOpDomainRhsPipeline());
702  CHKERR add_domain_ops_rhs_constrain(pipeline_mng->getOpDomainRhsPipeline());
703  CHKERR add_boundary_ops_rhs_mechanical(
704  pipeline_mng->getOpBoundaryRhsPipeline());
705 
706  CHKERR add_domain_base_ops(feAxiatorRhs->getOpPtrVector());
707  CHKERR add_domain_stress_ops(feAxiatorRhs->getOpPtrVector(),
708  commonPlasticDataPtr->mDPtr_Axiator);
709  CHKERR add_domain_ops_rhs_mechanical(feAxiatorRhs->getOpPtrVector());
710 
711  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule_deviator);
712  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule_deviator);
713 
714  CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule_bc);
715  CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule_bc);
716 
717  auto create_reaction_pipeline = [&](auto &pipeline) {
719 
720  if (reactionMarker) {
721 
722  if (SPACE_DIM == 2) {
723  auto det_ptr = boost::make_shared<VectorDouble>();
724  auto jac_ptr = boost::make_shared<MatrixDouble>();
725  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
726  pipeline.push_back(new OpCalculateHOJacForFace(jac_ptr));
727  pipeline.push_back(
728  new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
729  pipeline.push_back(new OpSetInvJacH1ForFace(inv_jac_ptr));
730  }
731 
732  pipeline.push_back(
734  "U", commonPlasticDataPtr->mGradPtr));
736  "EP", commonPlasticDataPtr->getPlasticStrainPtr()));
737 
738  if (is_large_strains) {
739 
740  if (commonPlasticDataPtr->mGradPtr != commonHenckyDataPtr->matGradPtr)
741  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
742  "Wrong pointer for grad");
743 
744  pipeline.push_back(
745  new OpCalculateEigenVals<SPACE_DIM>("U", commonHenckyDataPtr));
746  pipeline.push_back(
747  new OpCalculateLogC<SPACE_DIM>("U", commonHenckyDataPtr));
748  pipeline.push_back(
749  new OpCalculateLogC_dC<SPACE_DIM>("U", commonHenckyDataPtr));
750  pipeline.push_back(new OpCalculateHenckyPlasticStress<SPACE_DIM>(
751  "U", commonHenckyDataPtr, commonPlasticDataPtr->mDPtr));
752  pipeline.push_back(
753  new OpCalculatePiolaStress<SPACE_DIM>("U", commonHenckyDataPtr));
754 
755  } else {
756  pipeline.push_back(new OpSymmetrizeTensor<SPACE_DIM>(
757  "U", commonPlasticDataPtr->mGradPtr,
758  commonPlasticDataPtr->mStrainPtr));
759  pipeline.push_back(new OpPlasticStress("U", commonPlasticDataPtr,
760  commonPlasticDataPtr->mDPtr, 1));
761  }
762 
763  pipeline.push_back(new OpSetBc("U", false, reactionMarker));
764  // Calculate internal forece
765  if (is_large_strains) {
766  pipeline.push_back(new OpInternalForcePiola(
767  "U", commonHenckyDataPtr->getMatFirstPiolaStress()));
768  } else {
769  pipeline.push_back(
770  new OpInternalForceCauchy("U", commonPlasticDataPtr->mStressPtr));
771  }
772  pipeline.push_back(new OpUnSetBc("U"));
773  }
774 
776  };
777 
778  reactionFe = boost::make_shared<DomainEle>(mField);
779  reactionFe->getRuleHook = integration_rule_deviator;
780 
781  CHKERR create_reaction_pipeline(reactionFe->getOpPtrVector());
782 
784 }
785 //! [Push operators to pipeline]
786 
787 //! [Solve]
790 
791  Simple *simple = mField.getInterface<Simple>();
792  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
793  ISManager *is_manager = mField.getInterface<ISManager>();
794 
795  auto set_section_monitor = [&](auto solver) {
797  SNES snes;
798  CHKERR TSGetSNES(solver, &snes);
799  PetscViewerAndFormat *vf;
800  CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
801  PETSC_VIEWER_DEFAULT, &vf);
802  CHKERR SNESMonitorSet(
803  snes,
804  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
805  vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
807  };
808 
809  auto create_post_process_element = [&]() {
811  postProcFe = boost::make_shared<PostProcEle>(mField);
812  postProcFe->generateReferenceElementMesh();
813  if (SPACE_DIM == 2) {
814  auto det_ptr = boost::make_shared<VectorDouble>();
815  auto jac_ptr = boost::make_shared<MatrixDouble>();
816  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
817  postProcFe->getOpPtrVector().push_back(
818  new OpCalculateHOJacForFace(jac_ptr));
819  postProcFe->getOpPtrVector().push_back(
820  new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
821  postProcFe->getOpPtrVector().push_back(new OpSetInvJacH1ForFace(inv_jac_ptr));
822  }
823 
824  postProcFe->getOpPtrVector().push_back(
826  "U", commonPlasticDataPtr->mGradPtr));
827  postProcFe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
828  "TAU", commonPlasticDataPtr->getPlasticTauPtr()));
829  postProcFe->getOpPtrVector().push_back(
831  "EP", commonPlasticDataPtr->getPlasticStrainPtr()));
832 
833  if (is_large_strains) {
834 
835  if (commonPlasticDataPtr->mGradPtr != commonHenckyDataPtr->matGradPtr)
836  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
837  "Wrong pointer for grad");
838 
839  postProcFe->getOpPtrVector().push_back(
840  new OpCalculateEigenVals<SPACE_DIM>("U", commonHenckyDataPtr));
841  postProcFe->getOpPtrVector().push_back(
842  new OpCalculateLogC<SPACE_DIM>("U", commonHenckyDataPtr));
843  postProcFe->getOpPtrVector().push_back(
844  new OpCalculateLogC_dC<SPACE_DIM>("U", commonHenckyDataPtr));
845  postProcFe->getOpPtrVector().push_back(
847  "U", commonHenckyDataPtr, commonPlasticDataPtr->mDPtr, scale));
848  postProcFe->getOpPtrVector().push_back(
849  new OpCalculatePiolaStress<SPACE_DIM>("U", commonHenckyDataPtr));
850  postProcFe->getOpPtrVector().push_back(new OpPostProcHencky<SPACE_DIM>(
851  "U", postProcFe->postProcMesh, postProcFe->mapGaussPts,
852  commonHenckyDataPtr));
853 
854  } else {
855  postProcFe->getOpPtrVector().push_back(
856  new OpSymmetrizeTensor<SPACE_DIM>("U", commonPlasticDataPtr->mGradPtr,
857  commonPlasticDataPtr->mStrainPtr));
858  postProcFe->getOpPtrVector().push_back(new OpPlasticStress(
859  "U", commonPlasticDataPtr, commonPlasticDataPtr->mDPtr, scale));
860  postProcFe->getOpPtrVector().push_back(
862  "U", postProcFe->postProcMesh, postProcFe->mapGaussPts,
863  commonPlasticDataPtr->mStrainPtr,
864  commonPlasticDataPtr->mStressPtr));
865  }
866 
867  postProcFe->getOpPtrVector().push_back(
868  new OpCalculatePlasticSurface("U", commonPlasticDataPtr));
869  postProcFe->getOpPtrVector().push_back(
870  new OpPostProcPlastic("U", postProcFe->postProcMesh,
871  postProcFe->mapGaussPts, commonPlasticDataPtr));
872  postProcFe->addFieldValuesPostProc("U");
874  };
875 
876  auto scatter_create = [&](auto D, auto coeff) {
878  CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
879  ROW, "U", coeff, coeff, is);
880  int loc_size;
881  CHKERR ISGetLocalSize(is, &loc_size);
882  Vec v;
883  CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
884  VecScatter scatter;
885  CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
886  return std::make_tuple(SmartPetscObj<Vec>(v),
887  SmartPetscObj<VecScatter>(scatter));
888  };
889 
890  auto set_time_monitor = [&](auto dm, auto solver) {
892  boost::shared_ptr<Monitor> monitor_ptr(new Monitor(
893  dm, postProcFe, reactionFe, uXScatter, uYScatter, uZScatter));
894  boost::shared_ptr<ForcesAndSourcesCore> null;
895  CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
896  monitor_ptr, null, null);
898  };
899 
900  auto set_fieldsplit_preconditioner = [&](auto solver) {
902 
903  SNES snes;
904  CHKERR TSGetSNES(solver, &snes);
905  KSP ksp;
906  CHKERR SNESGetKSP(snes, &ksp);
907  PC pc;
908  CHKERR KSPGetPC(ksp, &pc);
909  PetscBool is_pcfs = PETSC_FALSE;
910  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
911 
912  // Setup fieldsplit (block) solver - optional: yes/no
913  if (is_pcfs == PETSC_TRUE) {
914 
915  auto bc_mng = mField.getInterface<BcManager>();
916  auto name_prb = simple->getProblemName();
917  auto is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_X", "U", 0, 0);
918  is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_Y", "U", 1, 1, is_all_bc);
919  is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_Z", "U", 2, 2, is_all_bc);
920  is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_ALL", "U", 0, 2, is_all_bc);
921 
922  int is_all_bc_size;
923  CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
924  MOFEM_LOG("EXAMPLE", Sev::inform)
925  << "Field split block size " << is_all_bc_size;
926 
927  CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
928  is_all_bc); // boundary block
929  }
930 
932  };
933 
934  auto dm = simple->getDM();
935  auto D = smartCreateDMVector(dm);
936 
937  boost::shared_ptr<FEMethod> null;
938  CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), feAxiatorLhs,
939  null, null);
940  CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(), feAxiatorRhs,
941  null, null);
942 
943  CHKERR create_post_process_element();
944  uXScatter = scatter_create(D, 0);
945  uYScatter = scatter_create(D, 1);
946  if (SPACE_DIM == 3)
947  uZScatter = scatter_create(D, 2);
948 
949  auto solver = pipeline_mng->createTSIM();
950 
951  CHKERR TSSetSolution(solver, D);
952  CHKERR set_section_monitor(solver);
953  CHKERR set_time_monitor(dm, solver);
954  CHKERR TSSetSolution(solver, D);
955  CHKERR TSSetFromOptions(solver);
956  CHKERR set_fieldsplit_preconditioner(solver);
957  CHKERR TSSetUp(solver);
958  CHKERR TSSolve(solver, NULL);
959 
960  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
961  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
962  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
963 
965 }
966 //! [Solve]
967 
968 static char help[] = "...\n\n";
969 
970 int main(int argc, char *argv[]) {
971 
972  // Initialisation of MoFEM/PETSc and MOAB data structures
973  const char param_file[] = "param_file.petsc";
974  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
975 
976  // Add logging channel for example
977  auto core_log = logging::core::get();
978  core_log->add_sink(
980  LogManager::setLog("EXAMPLE");
981  MOFEM_LOG_TAG("EXAMPLE", "example");
982 
983  try {
984 
985  //! [Register MoFEM discrete manager in PETSc]
986  DMType dm_name = "DMMOFEM";
987  CHKERR DMRegister_MoFEM(dm_name);
988  //! [Register MoFEM discrete manager in PETSc
989 
990  //! [Create MoAB]
991  moab::Core mb_instance; ///< mesh database
992  moab::Interface &moab = mb_instance; ///< mesh database interface
993  //! [Create MoAB]
994 
995  //! [Create MoFEM]
996  MoFEM::Core core(moab); ///< finite element database
997  MoFEM::Interface &m_field = core; ///< finite element database insterface
998  //! [Create MoFEM]
999 
1000  //! [Load mesh]
1001  Simple *simple = m_field.getInterface<Simple>();
1002  CHKERR simple->getOptions();
1003  CHKERR simple->loadFile();
1004  //! [Load mesh]
1005 
1006  //! [Example]
1007  Example ex(m_field);
1008  CHKERR ex.runProblem();
1009  //! [Example]
1010  }
1011  CATCH_ERRORS;
1012 
1014 }
std::string param_file
constexpr double third
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
Kronecker Delta class symmetric.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
constexpr double shear_modulus_G
constexpr double bulk_modulus_K
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, SPACE_DIM > OpBodyForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
@ ROW
Definition: definitions.h:136
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
FieldApproximationBase
approximation base
Definition: definitions.h:71
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ H1
continuous field
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ BLOCKSET
Definition: definitions.h:161
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
@ MOFEM_INVALID_DATA
Definition: definitions.h:49
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
constexpr auto t_kd
void temp(int x, int y=10)
Definition: simple.cpp:4
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:956
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
Definition: DMMMoFEM.cpp:755
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:478
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
Definition: DMMMoFEM.cpp:808
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:364
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:309
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:340
#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.
FTensor::Index< 'i', SPACE_DIM > i
double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition: DMMMoFEM.cpp:991
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
CoreTmp< 0 > Core
Definition: Core.hpp:1095
FaceElementForcesAndSourcesCore BoundaryEle
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1949
const double D
diffusivity
double A
double young_modulus
Definition: plastic.cpp:105
int main(int argc, char *argv[])
Definition: plastic.cpp:970
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: plastic.cpp:88
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: plastic.cpp:86
double Qinf
Definition: plastic.cpp:112
static char help[]
[Solve]
Definition: plastic.cpp:968
double rho
Definition: plastic.cpp:107
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: plastic.cpp:74
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:24
constexpr int SPACE_DIM
Definition: plastic.cpp:51
double visH
Definition: plastic.cpp:110
double poisson_ratio
Definition: plastic.cpp:106
double scale
Definition: plastic.cpp:103
long double hardening(long double tau, double temp)
Definition: plastic.cpp:116
double H
Definition: plastic.cpp:109
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
Definition: plastic.cpp:95
DataForcesAndSourcesCore::EntData EntData
Definition: plastic.cpp:55
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Body force]
Definition: plastic.cpp:72
long double hardening_dtau(long double tau, double temp)
Definition: plastic.cpp:120
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
Definition: plastic.cpp:97
int order
Definition: plastic.cpp:114
double b_iso
Definition: plastic.cpp:113
PetscBool is_large_strains
Definition: plastic.cpp:101
double sigmaY
Definition: plastic.cpp:108
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: plastic.cpp:93
double cn
Definition: plastic.cpp:111
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, SPACE_DIM > OpBodyForce
[Body force]
Definition: plastic.cpp:67
constexpr EntityType boundary_ent
Definition: plastic.cpp:54
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:41
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
PipelineManager::FaceEle DomainEle
[Example]
boost::shared_ptr< DomainEle > feAxiatorRhs
Definition: plastic.cpp:149
boost::shared_ptr< DomainEle > reactionFe
Definition: plastic.cpp:148
boost::shared_ptr< std::vector< unsigned char > > reactionMarker
Definition: plastic.cpp:157
boost::shared_ptr< PostProcEle > postProcFe
Definition: plastic.cpp:147
boost::shared_ptr< DomainEle > feAxiatorLhs
Definition: plastic.cpp:150
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
Definition: plastic.cpp:788
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: plastic.cpp:153
MoFEMErrorCode createCommonData()
std::vector< FTensor::Tensor1< double, 3 > > bodyForces
Definition: plastic.cpp:159
Example(MoFEM::Interface &m_field)
Definition: plastic.cpp:132
MoFEMErrorCode OPs()
[Boundary condition]
Definition: plastic.cpp:387
MoFEMErrorCode runProblem()
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: plastic.cpp:154
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Definition: plastic.cpp:156
MoFEMErrorCode setupProblem()
MoFEMErrorCode bC()
[Create common data]
Definition: plastic.cpp:319
boost::shared_ptr< PlasticOps::CommonData > commonPlasticDataPtr
Definition: plastic.cpp:145
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: plastic.cpp:152
boost::shared_ptr< HenckyOps::CommonData > commonHenckyDataPtr
Definition: plastic.cpp:146
Definition: HenckyOps.hpp:512
Simple interface for fast problem set-up.
Definition: BcManager.hpp:27
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
Definition: BcManager.hpp:131
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Data on single entity (This is passed as argument to DataOperator::doWork)
Deprecated interface functions.
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:33
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:279
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:323
Get value at integration points for scalar field.
Calculate symmetric tensor field rates ant integratio pts.
Calculate symmetric tensor field values at integration pts.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Scale base functions by inverses of measure of element.
Set indices on entities on finite element.
PipelineManager interface.
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:33
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
Definition: PlasticOps.hpp:279
Postprocess on face.
Post processing.
[Class definition]