v0.13.0
thermo_plastic.cpp
Go to the documentation of this file.
1 /**
2  * \file thermo_plastic.cpp
3  * \example thermo_plastic.cpp
4  *
5  * Thermo plasticity
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 // Thermal operators
102 /**
103  * @brief Integrate Lhs base of flux (1/k) base of flux (FLUX x FLUX)
104  *
105  */
108 
109 /**
110  * @brief Integrate Lhs div of base of flux time base of temperature (FLUX x T)
111  * and transpose of it, i.e. (T x FLAX)
112  *
113  */
115  GAUSS>::OpMixDivTimesScalar<SPACE_DIM>;
116 
117 /**
118  * @brief Integrate Lhs base of temerature times (heat capacity) times base of
119  * temperature (T x T)
120  *
121  */
124 
125 /**
126  * @brief Integrating Rhs flux base (1/k) flux (FLUX)
127  */
129  GAUSS>::OpBaseTimesVector<3, 3, 1>;
130 
131 /**
132  * @brief Integrate Rhs div flux base times temperature (T)
133  *
134  */
136  GAUSS>::OpMixDivTimesU<3, 1, 2>;
137 
138 /**
139  * @brief Integrate Rhs base of temerature time heat capacity times heat rate
140  * (T)
141  *
142  */
145 
146 /**
147  * @brief Integrate Rhs base of temerature times divergenc of flux (T)
148  *
149  */
151 
153  GAUSS>::OpSource<1, 1>;
156 
157 double scale = 1.;
158 
159 double young_modulus = 206913;
160 double poisson_ratio = 0.29;
161 double rho = 0;
162 double sigmaY = 450;
163 double H = 129;
164 double visH = 0;
165 double cn = 1;
166 double Qinf = 0; // 265;
167 double b_iso = 16.93;
169  16.2; // Force / (time temerature ) or Power /
170  // (length temperature) // Time unit is hour. force unit kN
172  5961.6; // length^2/(time^2 temerature) // length is millimeter time is hour
173 double omega_0 = 2e-3;
174 double omega_h = 2e-3;
175 double omega_inf = 0;
177 int order = 2;
178 
180 double amplitude_cycle = 0.5;
181 double amplitude_shift = 0.5;
182 double phase_shift = 0.8;
183 
184 inline long double hardening(long double tau, double temp) {
185  return H * tau * (1 - omega_h * temp) +
186  Qinf * (1. - std::exp(-b_iso * tau)) * (1 - omega_inf * temp) +
187  sigmaY * (1 - omega_0 * temp);
188 }
189 
190 inline long double hardening_dtau(long double tau, double temp) {
191  return H * (1 - omega_h * temp) +
192  Qinf * b_iso * std::exp(-b_iso * tau) * (1 - omega_inf * temp);
193 }
194 
195 inline long double hardening_dtemp(long double tau, double temp) {
196  return -H * tau * omega_h - Qinf * (1. - std::exp(-b_iso * tau)) * omega_inf -
197  sigmaY * omega_0;
198 }
199 
200 #include <HenckyOps.hpp>
201 #include <PlasticOps.hpp>
202 #include <OpPostProcElastic.hpp>
203 using namespace PlasticOps;
204 #include <PlasticThermalOps.hpp>
205 
206 using namespace HenckyOps;
207 struct Example {
208 
209  Example(MoFEM::Interface &m_field) : mField(m_field) {}
210 
212 
213 private:
214  MoFEM::Interface &mField;
215 
221 
222  boost::shared_ptr<PlasticThermalOps::CommonData> commonPlasticDataPtr;
223  boost::shared_ptr<HenckyOps::CommonData> commonHenckyDataPtr;
224  boost::shared_ptr<PostProcEle> postProcFe;
225  boost::shared_ptr<DomainEle> reactionFe;
226  boost::shared_ptr<DomainEle> feAxiatorRhs;
227  boost::shared_ptr<DomainEle> feAxiatorLhs;
228  boost::shared_ptr<DomainEle> feThermalRhs;
229  boost::shared_ptr<DomainEle> feThermalLhs;
230 
231  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
232  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
233  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
234 
235  boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
236  boost::shared_ptr<std::vector<unsigned char>> reactionMarker;
237 
238  std::vector<FTensor::Tensor1<double, 3>> bodyForces;
239 
240  struct BcTempFun {
241  BcTempFun(double v, FEMethod &fe) : valTemp(v), fE(fe) {}
242  double operator()(const double, const double, const double) {
243  return -valTemp; // * fE.ts_t;
244  }
245 
246  private:
247  double valTemp;
249  };
250  std::vector<BcTempFun> bcTemperatureFunctions;
251 };
252 
253 //! [Run problem]
256  CHKERR setupProblem();
257  CHKERR createCommonData();
258  CHKERR bC();
259  CHKERR OPs();
260  CHKERR tsSolve();
262 }
263 //! [Run problem]
264 
265 //! [Set up problem]
268  Simple *simple = mField.getInterface<Simple>();
269  // Add field
271  constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
272  // Mechanical fields
273  CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
274  CHKERR simple->addDomainField("TAU", L2, base, 1);
275  CHKERR simple->addDomainField("EP", L2, base, size_symm);
276  CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
277  // Temerature
278  const auto flux_space = (SPACE_DIM == 2) ? HCURL : HDIV;
279  CHKERR simple->addDomainField("T", L2, AINSWORTH_LEGENDRE_BASE, 1);
280  CHKERR simple->addDomainField("FLUX", flux_space, DEMKOWICZ_JACOBI_BASE, 1);
281  CHKERR simple->addBoundaryField("FLUX", flux_space, DEMKOWICZ_JACOBI_BASE, 1);
282 
283  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
284  CHKERR simple->setFieldOrder("U", order);
285  CHKERR simple->setFieldOrder("TAU", order - 1);
286  CHKERR simple->setFieldOrder("EP", order - 1);
287  CHKERR simple->setFieldOrder("FLUX", order);
288  CHKERR simple->setFieldOrder("T", order - 1);
289  CHKERR simple->setUp();
291 }
292 //! [Set up problem]
293 
294 //! [Create common data]
297 
298  auto get_command_line_parameters = [&]() {
300  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-scale", &scale, PETSC_NULL);
301  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-rho", &rho, PETSC_NULL);
302  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
303  &young_modulus, PETSC_NULL);
304  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
305  &poisson_ratio, PETSC_NULL);
306  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening", &H, PETSC_NULL);
307  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening_viscous", &visH,
308  PETSC_NULL);
309  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-yield_stress", &sigmaY,
310  PETSC_NULL);
311  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn", &cn, PETSC_NULL);
312  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-Qinf", &Qinf, PETSC_NULL);
313  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-b_iso", &b_iso, PETSC_NULL);
314 
315  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-capacity", &heat_capacity,
316  PETSC_NULL);
317  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-conductivity",
318  &heat_conductivity, PETSC_NULL);
319  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-omega_0", &omega_0,
320  PETSC_NULL);
321  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-omega_h", &omega_0,
322  PETSC_NULL);
323  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-omega_inf", &omega_inf,
324  PETSC_NULL);
325  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-fraction_of_dissipation",
326  &fraction_of_dissipation, PETSC_NULL);
327  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-number_of_cycles_in_total_time",
328  &number_of_cycles_in_total_time, PETSC_NULL);
329  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-amplitude_cycle",
330  &amplitude_cycle, PETSC_NULL);
331  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-phase_shift", &phase_shift,
332  PETSC_NULL);
333  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-amplitude_shift",
334  &amplitude_shift, PETSC_NULL);
335  MOFEM_LOG("EXAMPLE", Sev::inform) << "Young modulus " << young_modulus;
336  MOFEM_LOG("EXAMPLE", Sev::inform) << "Poisson ratio " << poisson_ratio;
337  MOFEM_LOG("EXAMPLE", Sev::inform) << "Yield stress " << sigmaY;
338  MOFEM_LOG("EXAMPLE", Sev::inform) << "Hardening " << H;
339  MOFEM_LOG("EXAMPLE", Sev::inform) << "Viscous hardening " << visH;
340  MOFEM_LOG("EXAMPLE", Sev::inform) << "Saturation yield stress " << Qinf;
341  MOFEM_LOG("EXAMPLE", Sev::inform) << "Saturation exponent " << b_iso;
342  MOFEM_LOG("EXAMPLE", Sev::inform) << "cn " << cn;
343  MOFEM_LOG("EXAMPLE", Sev::inform)
344  << "fraction_of_dissipation " << fraction_of_dissipation;
345 
346  PetscBool is_scale = PETSC_TRUE;
347  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_scale", &is_scale,
348  PETSC_NULL);
349  if (is_scale) {
351  young_modulus *= scale;
352  rho *= scale;
353  sigmaY *= scale;
354  H *= scale;
355  Qinf *= scale;
356  visH *= scale;
357  heat_capacity *= scale;
358 
359  MOFEM_LOG("EXAMPLE", Sev::inform)
360  << "Scaled Young modulus " << young_modulus;
361  MOFEM_LOG("EXAMPLE", Sev::inform)
362  << "Scaled Poisson ratio " << poisson_ratio;
363  MOFEM_LOG("EXAMPLE", Sev::inform) << "Scaled Yield stress " << sigmaY;
364  MOFEM_LOG("EXAMPLE", Sev::inform) << "Scaled Hardening " << H;
365  MOFEM_LOG("EXAMPLE", Sev::inform) << "Scaled Viscous hardening " << visH;
366  MOFEM_LOG("EXAMPLE", Sev::inform)
367  << "Scaled Saturation yield stress " << Qinf;
368  MOFEM_LOG("EXAMPLE", Sev::inform)
369  << "Scaled heat capacity " << heat_capacity;
370  }
371 
373  };
374 
375  auto set_matrial_stiffness = [&]() {
382  const double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
383  const double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
384 
385  // Plane stress or when 1, plane strain or 3d
386  const double A = (SPACE_DIM == 2)
387  ? 2 * shear_modulus_G /
388  (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
389  : 1;
390 
391  auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
392  *commonPlasticDataPtr->mDPtr);
393  auto t_D_axiator = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
394  *commonPlasticDataPtr->mDPtr_Axiator);
395  auto t_D_deviator = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
396  *commonPlasticDataPtr->mDPtr_Deviator);
397 
398  constexpr double third = boost::math::constants::third<double>();
399  t_D_axiator(i, j, k, l) = A *
400  (bulk_modulus_K - (2. / 3.) * shear_modulus_G) *
401  t_kd(i, j) * t_kd(k, l);
402  t_D_deviator(i, j, k, l) =
403  2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
404  t_D(i, j, k, l) = t_D_axiator(i, j, k, l) + t_D_deviator(i, j, k, l);
405 
407  };
408 
409  commonPlasticDataPtr = boost::make_shared<PlasticThermalOps::CommonData>();
410  constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
411  commonPlasticDataPtr->mDPtr = boost::make_shared<MatrixDouble>();
412  commonPlasticDataPtr->mDPtr->resize(size_symm * size_symm, 1);
413  commonPlasticDataPtr->mDPtr_Axiator = boost::make_shared<MatrixDouble>();
414  commonPlasticDataPtr->mDPtr_Axiator->resize(size_symm * size_symm, 1);
415  commonPlasticDataPtr->mDPtr_Deviator = boost::make_shared<MatrixDouble>();
416  commonPlasticDataPtr->mDPtr_Deviator->resize(size_symm * size_symm, 1);
417 
418  commonPlasticDataPtr->mGradPtr = boost::make_shared<MatrixDouble>();
419  commonPlasticDataPtr->mStrainPtr = boost::make_shared<MatrixDouble>();
420  commonPlasticDataPtr->mStressPtr = boost::make_shared<MatrixDouble>();
421 
422  CHKERR get_command_line_parameters();
423  CHKERR set_matrial_stiffness();
424 
426 }
427 //! [Create common data]
428 
429 //! [Boundary condition]
432 
433  auto simple = mField.getInterface<Simple>();
434  auto bc_mng = mField.getInterface<BcManager>();
435  auto prb_mng = mField.getInterface<ProblemsManager>();
436 
437  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
438  "U", 0, 0);
439  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
440  "U", 1, 1);
441  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
442  "U", 2, 2);
443  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
444  "REMOVE_ALL", "U", 0, 3);
445  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
446  "ZERO_FLUX", "FLUX", 0, 1);
447 
448  PetscBool zero_fix_skin_flux = PETSC_TRUE;
449  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-fix_skin_flux",
450  &zero_fix_skin_flux, PETSC_NULL);
451  if (zero_fix_skin_flux) {
452  Range faces;
453  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, faces);
454  Skinner skin(&mField.get_moab());
455  Range skin_ents;
456  CHKERR skin.find_skin(0, faces, false, skin_ents);
457  ParallelComm *pcomm =
458  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
459  if (pcomm == NULL)
460  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
461  "Communicator not created");
462  CHKERR pcomm->filter_pstatus(skin_ents,
463  PSTATUS_SHARED | PSTATUS_MULTISHARED,
464  PSTATUS_NOT, -1, &skin_ents);
465  CHKERR prb_mng->removeDofsOnEntities(simple->getProblemName(), "FLUX",
466  skin_ents, 0, 1);
467  }
468 
469  CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_X", "U",
470  0, 0);
471  CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_Y", "U",
472  1, 1);
473  CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_Z", "U",
474  2, 2);
475  CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
476  "U", 0, 3);
477 
478  auto &bc_map = bc_mng->getBcMapByBlockName();
479  boundaryMarker = bc_mng->getMergedBlocksMarker(vector<string>{"FIX_"});
480 
481  CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "REACTION",
482  "U", 0, 3);
483  for (auto b : bc_map) {
484  MOFEM_LOG("EXAMPLE", Sev::verbose) << "Marker " << b.first;
485  }
486 
487  // OK. We have problem with GMesh, it adding empty characters at the end of
488  // block. So first block is search by regexp. popMarkDOFsOnEntities should
489  // work with regexp.
490  std::string reaction_block_set;
491  for (auto bc : bc_map) {
492  if (bc_mng->checkBlock(bc, "REACTION")) {
493  reaction_block_set = bc.first;
494  break;
495  }
496  }
497 
498  if (auto bc = bc_mng->popMarkDOFsOnEntities(reaction_block_set)) {
499  reactionMarker = bc->getBcMarkersPtr();
500 
501  // Only take reaction from nodes
502  Range nodes;
503  CHKERR mField.get_moab().get_entities_by_type(0, MBVERTEX, nodes, true);
504  CHKERR prb_mng->markDofs(simple->getProblemName(), ROW,
505  ProblemsManager::MarkOP::AND, nodes,
506  *reactionMarker);
507 
508  } else {
509  MOFEM_LOG("EXAMPLE", Sev::warning) << "REACTION blockset does not exist";
510  }
511 
512  if (!reactionMarker) {
513  MOFEM_LOG("EXAMPLE", Sev::warning) << "REACTION blockset does not exist";
514  }
515 
517 }
518 //! [Boundary condition]
519 
520 //! [Push operators to pipeline]
523  auto pipeline_mng = mField.getInterface<PipelineManager>();
524  auto simple = mField.getInterface<Simple>();
525  auto bc_mng = mField.getInterface<BcManager>();
526 
527  auto integration_rule_deviator = [](int o_row, int o_col, int approx_order) {
528  return 2 * (approx_order - 1);
529  };
530  auto integration_rule_bc = [](int, int, int approx_order) {
531  return 2 * approx_order;
532  };
533 
534  feAxiatorRhs = boost::make_shared<DomainEle>(mField);
535  feAxiatorLhs = boost::make_shared<DomainEle>(mField);
536  auto integration_rule_axiator = [](int, int, int approx_order) {
537  return 2 * (approx_order - 1);
538  };
539  feAxiatorRhs->getRuleHook = integration_rule_axiator;
540  feAxiatorLhs->getRuleHook = integration_rule_axiator;
541 
542  feThermalRhs = boost::make_shared<DomainEle>(mField);
543  feThermalLhs = boost::make_shared<DomainEle>(mField);
544  auto integration_rule_thermal = [](int, int, int approx_order) {
545  return 2 * approx_order;
546  };
547  feThermalRhs->getRuleHook = integration_rule_thermal;
548  feThermalLhs->getRuleHook = integration_rule_thermal;
549 
550  auto add_domain_base_ops = [&](auto &pipeline) {
552 
553  if (SPACE_DIM == 2) {
554  auto det_ptr = boost::make_shared<VectorDouble>();
555  auto jac_ptr = boost::make_shared<MatrixDouble>();
556  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
557  pipeline.push_back(new OpCalculateHOJacForFace(jac_ptr));
558  pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
559  pipeline.push_back(new OpSetInvJacH1ForFace(inv_jac_ptr));
560  pipeline.push_back(new OpMakeHdivFromHcurl());
561  pipeline.push_back(new OpSetContravariantPiolaTransformOnFace2D(jac_ptr));
562  pipeline.push_back(new OpSetInvJacHcurlFace(inv_jac_ptr));
563  pipeline.push_back(new OpSetInvJacL2ForFace(inv_jac_ptr));
564  }
565 
566  pipeline.push_back(new OpCalculateScalarFieldValuesDot(
567  "TAU", commonPlasticDataPtr->getPlasticTauDotPtr()));
569  "EP", commonPlasticDataPtr->getPlasticStrainPtr()));
571  "EP", commonPlasticDataPtr->getPlasticStrainDotPtr()));
573  "U", commonPlasticDataPtr->mGradPtr));
574  pipeline.push_back(new OpCalculateScalarFieldValues(
575  "TAU", commonPlasticDataPtr->getPlasticTauPtr()));
576 
577  pipeline.push_back(new OpCalculateScalarFieldValues(
578  "T", commonPlasticDataPtr->getTempValPtr()));
579  pipeline.push_back(new OpCalculateScalarFieldValuesDot(
580  "T", commonPlasticDataPtr->getTempValDotPtr()));
581  pipeline.push_back(new OpCalculateHdivVectorDivergence<3, SPACE_DIM>(
582  "FLUX", commonPlasticDataPtr->getTempDivFluxPtr()));
583  pipeline.push_back(new OpCalculateHVecVectorField<3>(
584  "FLUX", commonPlasticDataPtr->getTempFluxValPtr()));
585 
587  };
588 
589  auto add_domain_stress_ops = [&](auto &pipeline, auto m_D_ptr) {
591 
592  pipeline.push_back(new OpSymmetrizeTensor<SPACE_DIM>(
593  "U", commonPlasticDataPtr->mGradPtr, commonPlasticDataPtr->mStrainPtr));
594  pipeline.push_back(
595  new OpPlasticStress("U", commonPlasticDataPtr, m_D_ptr, 1));
596 
597  if (m_D_ptr != commonPlasticDataPtr->mDPtr_Axiator)
598  pipeline.push_back(
599  new OpCalculatePlasticSurface("U", commonPlasticDataPtr));
600 
602  };
603 
604  auto add_domain_ops_lhs_mechanical = [&](auto &pipeline, auto m_D_ptr) {
606  pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
607 
608  pipeline.push_back(new OpKCauchy("U", "U", m_D_ptr));
609  pipeline.push_back(new OpCalculatePlasticInternalForceLhs_dEP(
610  "U", "EP", commonPlasticDataPtr, m_D_ptr));
611 
612  pipeline.push_back(new OpUnSetBc("U"));
614  };
615 
616  auto add_domain_ops_rhs_mechanical = [&](auto &pipeline) {
618  pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
619 
621  const std::string block_name = "BODY_FORCE";
622  if (it->getName().compare(0, block_name.size(), block_name) == 0) {
623  std::vector<double> attr;
624  CHKERR it->getAttributes(attr);
625  if (attr.size() == 3) {
626  bodyForces.push_back(
627  FTensor::Tensor1<double, 3>{attr[0], attr[1], attr[2]});
628  } else {
629  SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
630  "Should be three atributes in BODYFORCE blockset, but is %d",
631  attr.size());
632  }
633  }
634  }
635 
636  auto get_body_force = [this](const double, const double, const double) {
637  auto *pipeline_mng = mField.getInterface<PipelineManager>();
640  t_source(i) = 0;
641  auto fe_domain_rhs = pipeline_mng->getDomainRhsFE();
642  const auto time = fe_domain_rhs->ts_t;
643  // hardcoded gravity load
644  for (auto &t_b : bodyForces)
645  t_source(i) += (scale * t_b(i)) * time;
646  return t_source;
647  };
648 
649  pipeline.push_back(new OpBodyForce("U", get_body_force));
650 
651  // Calculate internal forece
652  pipeline.push_back(
653  new OpInternalForceCauchy("U", commonPlasticDataPtr->mStressPtr));
654 
655  pipeline.push_back(new OpUnSetBc("U"));
657  };
658 
659  auto add_domain_ops_lhs_constrain = [&](auto &pipeline, auto m_D_ptr) {
661  pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
662 
663  pipeline.push_back(new OpCalculatePlasticFlowLhs_dU(
664  "EP", "U", commonPlasticDataPtr, m_D_ptr));
665  pipeline.push_back(new OpCalculateContrainsLhs_dU(
666  "TAU", "U", commonPlasticDataPtr, m_D_ptr));
667 
668  pipeline.push_back(new OpCalculatePlasticFlowLhs_dEP(
669  "EP", "EP", commonPlasticDataPtr, m_D_ptr));
670  pipeline.push_back(
671  new OpCalculatePlasticFlowLhs_dTAU("EP", "TAU", commonPlasticDataPtr));
672  pipeline.push_back(new OpCalculateContrainsLhs_dEP(
673  "TAU", "EP", commonPlasticDataPtr, m_D_ptr));
674  pipeline.push_back(
675  new OpCalculateContrainsLhs_dTAU("TAU", "TAU", commonPlasticDataPtr));
676 
677  pipeline.push_back(new PlasticThermalOps::OpCalculateContrainsLhs_dT(
678  "TAU", "T", commonPlasticDataPtr));
679  pipeline.push_back(new PlasticThermalOps::OpPlasticHeatProduction_dEP(
680  "T", "EP", commonPlasticDataPtr));
681 
682  pipeline.push_back(new OpUnSetBc("U"));
684  };
685 
686  auto add_domain_ops_rhs_constrain = [&](auto &pipeline) {
688  pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
689 
690  pipeline.push_back(
691  new OpCalculatePlasticFlowRhs("EP", commonPlasticDataPtr));
692  pipeline.push_back(
693  new OpCalculateContrainsRhs("TAU", commonPlasticDataPtr));
694 
695  pipeline.push_back(new PlasticThermalOps::OpPlasticHeatProduction(
696  "T", commonPlasticDataPtr));
697 
699  };
700 
701  auto add_boundary_ops_lhs_mechanical = [&](auto &pipeline) {
703  auto &bc_map = mField.getInterface<BcManager>()->getBcMapByBlockName();
704  for (auto bc : bc_map) {
705  if (bc_mng->checkBlock(bc, "FIX_")) {
706  pipeline.push_back(
707  new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
708  pipeline.push_back(new OpBoundaryMass(
709  "U", "U", [](double, double, double) { return 1.; },
710  bc.second->getBcEntsPtr()));
711  pipeline.push_back(new OpUnSetBc("U"));
712  }
713  }
715  };
716 
717  auto add_boundary_ops_rhs_mechanical = [&](auto &pipeline) {
719 
720  auto get_time = [&](double, double, double) {
721  auto *pipeline_mng = mField.getInterface<PipelineManager>();
722  auto &fe_domain_rhs = pipeline_mng->getBoundaryRhsFE();
723  return fe_domain_rhs->ts_t;
724  };
725 
726  auto get_time_scaled = [&](double, double, double) {
727  auto *pipeline_mng = mField.getInterface<PipelineManager>();
728  auto &fe_domain_rhs = pipeline_mng->getBoundaryRhsFE();
729  return fe_domain_rhs->ts_t * scale;
730  };
731 
732  auto get_minus_time = [&](double, double, double) {
733  auto *pipeline_mng = mField.getInterface<PipelineManager>();
734  auto &fe_domain_rhs = pipeline_mng->getBoundaryRhsFE();
735  return -fe_domain_rhs->ts_t;
736  };
737 
738  auto time_scaled = [&](double, double, double) {
739  auto *pipeline_mng = mField.getInterface<PipelineManager>();
740  auto &fe_domain_rhs = pipeline_mng->getBoundaryRhsFE();
742  return amplitude_cycle * sin((2 * fe_domain_rhs->ts_t - phase_shift) *
745  } else {
746  return fe_domain_rhs->ts_t;
747  }
748  };
749 
750  pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
751 
753  if (it->getName().compare(0, 5, "FORCE") == 0) {
754  Range force_edges;
755  std::vector<double> attr_vec;
756  CHKERR it->getMeshsetIdEntitiesByDimension(
757  mField.get_moab(), SPACE_DIM - 1, force_edges, true);
758  it->getAttributes(attr_vec);
759  if (attr_vec.size() < SPACE_DIM)
760  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
761  "Wrong size of boundary attributes vector. Set right block "
762  "size attributes.");
763  auto force_vec_ptr = boost::make_shared<MatrixDouble>(SPACE_DIM, 1);
764  std::copy(&attr_vec[0], &attr_vec[SPACE_DIM],
765  force_vec_ptr->data().begin());
766  pipeline.push_back(
767  new OpBoundaryVec("U", force_vec_ptr, get_time_scaled,
768  boost::make_shared<Range>(force_edges)));
769  }
770  }
771 
772  pipeline.push_back(new OpUnSetBc("U"));
773 
774  auto u_mat_ptr = boost::make_shared<MatrixDouble>();
775  pipeline.push_back(
776  new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_mat_ptr));
777 
778  for (auto &bc : mField.getInterface<BcManager>()->getBcMapByBlockName()) {
779  if (bc_mng->checkBlock(bc, "FIX_")) {
780  pipeline.push_back(
781  new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
782  auto attr_vec = boost::make_shared<MatrixDouble>(SPACE_DIM, 1);
783  attr_vec->clear();
784  if (bc.second->bcAttributes.size() < SPACE_DIM)
785  SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
786  "Wrong size of boundary attributes vector. Set right block "
787  "size attributes. Size of attributes %d",
788  bc.second->bcAttributes.size());
789  std::copy(&bc.second->bcAttributes[0],
790  &bc.second->bcAttributes[SPACE_DIM],
791  attr_vec->data().begin());
792 
793  pipeline.push_back(new OpBoundaryVec("U", attr_vec, time_scaled,
794  bc.second->getBcEntsPtr()));
795  pipeline.push_back(new OpBoundaryInternal(
796  "U", u_mat_ptr, [](double, double, double) { return 1.; },
797  bc.second->getBcEntsPtr()));
798 
799  pipeline.push_back(new OpUnSetBc("U"));
800  }
801  }
802 
804  };
805 
806  auto add_domain_ops_lhs_thermal = [&](auto &pipeline, auto &fe) {
808  auto resistance = [](const double, const double, const double) {
809  return (1. / heat_conductivity);
810  };
811  auto capacity = [&](const double, const double, const double) {
812  return -heat_capacity * fe.ts_a;
813  };
814  auto unity = []() { return 1; };
815  pipeline.push_back(new OpHdivHdiv("FLUX", "FLUX", resistance));
816  pipeline.push_back(new OpHdivT("FLUX", "T", unity, true));
817  pipeline.push_back(new OpCapacity("T", "T", capacity));
819  };
820 
821  auto add_domain_ops_rhs_thermal = [&](auto &pipeline) {
823  auto resistance = [](const double, const double, const double) {
824  return (1. / heat_conductivity);
825  };
826  auto capacity = [&](const double, const double, const double) {
827  return -heat_capacity;
828  };
829  auto unity = [](const double, const double, const double) { return 1; };
830  pipeline.push_back(new OpHdivFlux(
831  "FLUX", commonPlasticDataPtr->getTempFluxValPtr(), resistance));
832  pipeline.push_back(
833  new OpHDivTemp("FLUX", commonPlasticDataPtr->getTempValPtr(), unity));
834  pipeline.push_back(new OpBaseDivFlux(
835  "T", commonPlasticDataPtr->getTempDivFluxPtr(), unity));
836  // auto source = [&](const double x, const double y, const double z) {
837  // return 1;
838  // };
839  // pipeline.push_back(new OpHeatSource("T", source));
840 
841  pipeline.push_back(new OpBaseDotT(
842  "T", commonPlasticDataPtr->getTempValDotPtr(), capacity));
844  };
845 
846  auto add_boundary_ops_rhs_thermal = [&](auto &pipeline, auto &fe) {
848 
849  if (SPACE_DIM == 2) {
850  pipeline.push_back(new OpSetContravariantPiolaTransformOnEdge2D());
851  } else if (SPACE_DIM == 3) {
852  pipeline.push_back(new OpHOSetCovariantPiolaTransformOnFace3D(HDIV));
853  }
854 
856  const std::string block_name = "TEMPERATURE";
857  if (it->getName().compare(0, block_name.size(), block_name) == 0) {
858  Range temp_edges;
859  std::vector<double> attr_vec;
860  CHKERR it->getMeshsetIdEntitiesByDimension(
861  mField.get_moab(), SPACE_DIM - 1, temp_edges, true);
862  it->getAttributes(attr_vec);
863  if (attr_vec.size() != 1)
864  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
865  "Should be one attribute");
866  MOFEM_LOG("EXAMPLE", Sev::inform)
867  << "Set temerature " << attr_vec[0] << " on ents:\n"
868  << temp_edges;
869  bcTemperatureFunctions.push_back(BcTempFun(attr_vec[0], fe));
870  pipeline.push_back(
871  new OpTemperatureBC("FLUX", bcTemperatureFunctions.back(),
872  boost::make_shared<Range>(temp_edges)));
873  }
874  }
876  };
877 
878  // Mechanics
879  CHKERR add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
880  CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainLhsPipeline(),
881  commonPlasticDataPtr->mDPtr_Deviator);
882  CHKERR add_domain_ops_lhs_mechanical(pipeline_mng->getOpDomainLhsPipeline(),
883  commonPlasticDataPtr->mDPtr_Deviator);
884  CHKERR add_domain_ops_lhs_constrain(pipeline_mng->getOpDomainLhsPipeline(),
885  commonPlasticDataPtr->mDPtr_Deviator);
886  CHKERR add_boundary_ops_lhs_mechanical(
887  pipeline_mng->getOpBoundaryLhsPipeline());
888 
889  // Axiator
890  CHKERR add_domain_base_ops(feAxiatorLhs->getOpPtrVector());
891  CHKERR add_domain_ops_lhs_mechanical(feAxiatorLhs->getOpPtrVector(),
892  commonPlasticDataPtr->mDPtr_Axiator);
893  // Temperature
894  CHKERR add_domain_base_ops(feThermalLhs->getOpPtrVector());
895  CHKERR add_domain_ops_lhs_thermal(feThermalLhs->getOpPtrVector(),
896  *feThermalLhs);
897 
898  // Mechanics
899  CHKERR add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
900  CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainRhsPipeline(),
901  commonPlasticDataPtr->mDPtr_Deviator);
902  CHKERR add_domain_ops_rhs_mechanical(pipeline_mng->getOpDomainRhsPipeline());
903  CHKERR add_domain_ops_rhs_constrain(pipeline_mng->getOpDomainRhsPipeline());
904  CHKERR add_boundary_ops_rhs_mechanical(
905  pipeline_mng->getOpBoundaryRhsPipeline());
906 
907  // Axiator
908  CHKERR add_domain_base_ops(feAxiatorRhs->getOpPtrVector());
909  CHKERR add_domain_stress_ops(feAxiatorRhs->getOpPtrVector(),
910  commonPlasticDataPtr->mDPtr_Axiator);
911  CHKERR add_domain_ops_rhs_mechanical(feAxiatorRhs->getOpPtrVector());
912 
913  // Temperature
914  CHKERR add_domain_base_ops(feThermalRhs->getOpPtrVector());
915  CHKERR add_domain_ops_rhs_thermal(feThermalRhs->getOpPtrVector());
916  CHKERR
917  add_boundary_ops_rhs_thermal(pipeline_mng->getOpBoundaryRhsPipeline(),
918  *pipeline_mng->getBoundaryRhsFE());
919 
920  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule_deviator);
921  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule_deviator);
922 
923  CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule_bc);
924  CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule_bc);
925 
926  auto create_reaction_pipeline = [&](auto &pipeline) {
928 
929  if (reactionMarker) {
930 
931  if (SPACE_DIM == 2) {
932  auto det_ptr = boost::make_shared<VectorDouble>();
933  auto jac_ptr = boost::make_shared<MatrixDouble>();
934  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
935  pipeline.push_back(new OpCalculateHOJacForFace(jac_ptr));
936  pipeline.push_back(
937  new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
938  pipeline.push_back(new OpSetInvJacH1ForFace(inv_jac_ptr));
939  }
940 
941  pipeline.push_back(
943  "U", commonPlasticDataPtr->mGradPtr));
945  "EP", commonPlasticDataPtr->getPlasticStrainPtr()));
946  pipeline.push_back(
947  new OpSymmetrizeTensor<SPACE_DIM>("U", commonPlasticDataPtr->mGradPtr,
948  commonPlasticDataPtr->mStrainPtr));
949  pipeline.push_back(new OpPlasticStress("U", commonPlasticDataPtr,
950  commonPlasticDataPtr->mDPtr, 1));
951  pipeline.push_back(new OpSetBc("U", false, reactionMarker));
952 
953  // Calculate internal forece
954  pipeline.push_back(
955  new OpInternalForceCauchy("U", commonPlasticDataPtr->mStressPtr));
956 
957  pipeline.push_back(new OpUnSetBc("U"));
958  }
959 
961  };
962 
963  reactionFe = boost::make_shared<DomainEle>(mField);
964  reactionFe->getRuleHook = integration_rule_deviator;
965 
966  CHKERR create_reaction_pipeline(reactionFe->getOpPtrVector());
967 
969 }
970 //! [Push operators to pipeline]
971 
972 //! [Solve]
975 
976  Simple *simple = mField.getInterface<Simple>();
977  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
978  ISManager *is_manager = mField.getInterface<ISManager>();
979 
980  auto set_section_monitor = [&](auto solver) {
982  SNES snes;
983  CHKERR TSGetSNES(solver, &snes);
984  PetscViewerAndFormat *vf;
985  CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
986  PETSC_VIEWER_DEFAULT, &vf);
987  CHKERR SNESMonitorSet(
988  snes,
989  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
990  vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
992  };
993 
994  auto create_post_process_element = [&]() {
996  postProcFe = boost::make_shared<PostProcEle>(mField);
997  postProcFe->generateReferenceElementMesh();
998  if (SPACE_DIM == 2) {
999  auto det_ptr = boost::make_shared<VectorDouble>();
1000  auto jac_ptr = boost::make_shared<MatrixDouble>();
1001  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
1002  postProcFe->getOpPtrVector().push_back(
1003  new OpCalculateHOJacForFace(jac_ptr));
1004  postProcFe->getOpPtrVector().push_back(
1005  new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
1006  postProcFe->getOpPtrVector().push_back(
1007  new OpSetInvJacH1ForFace(inv_jac_ptr));
1008  }
1009 
1010  postProcFe->getOpPtrVector().push_back(
1012  "U", commonPlasticDataPtr->mGradPtr));
1013  postProcFe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
1014  "TAU", commonPlasticDataPtr->getPlasticTauPtr()));
1015  postProcFe->getOpPtrVector().push_back(
1017  "EP", commonPlasticDataPtr->getPlasticStrainPtr()));
1018  postProcFe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
1019  "T", commonPlasticDataPtr->getTempValPtr()));
1020 
1021  postProcFe->getOpPtrVector().push_back(new OpSymmetrizeTensor<SPACE_DIM>(
1022  "U", commonPlasticDataPtr->mGradPtr, commonPlasticDataPtr->mStrainPtr));
1023  postProcFe->getOpPtrVector().push_back(new OpPlasticStress(
1024  "U", commonPlasticDataPtr, commonPlasticDataPtr->mDPtr, scale));
1025  postProcFe->getOpPtrVector().push_back(
1027  "U", postProcFe->postProcMesh, postProcFe->mapGaussPts,
1028  commonPlasticDataPtr->mStrainPtr,
1029  commonPlasticDataPtr->mStressPtr));
1030 
1031  postProcFe->getOpPtrVector().push_back(
1032  new OpCalculatePlasticSurface("U", commonPlasticDataPtr));
1033  postProcFe->getOpPtrVector().push_back(
1034  new OpPostProcPlastic("U", postProcFe->postProcMesh,
1035  postProcFe->mapGaussPts, commonPlasticDataPtr));
1036  postProcFe->addFieldValuesPostProc("U");
1038  };
1039 
1040  auto scatter_create = [&](auto D, auto coeff) {
1041  SmartPetscObj<IS> is;
1042  CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
1043  ROW, "U", coeff, coeff, is);
1044  int loc_size;
1045  CHKERR ISGetLocalSize(is, &loc_size);
1046  Vec v;
1047  CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
1048  VecScatter scatter;
1049  CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
1050  return std::make_tuple(SmartPetscObj<Vec>(v),
1051  SmartPetscObj<VecScatter>(scatter));
1052  };
1053 
1054  auto set_time_monitor = [&](auto dm, auto solver) {
1056  boost::shared_ptr<Monitor> monitor_ptr(new Monitor(
1057  dm, postProcFe, reactionFe, uXScatter, uYScatter, uZScatter));
1058  boost::shared_ptr<ForcesAndSourcesCore> null;
1059  CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
1060  monitor_ptr, null, null);
1062  };
1063 
1064  auto dm = simple->getDM();
1065  auto D = smartCreateDMVector(dm);
1066 
1067  boost::shared_ptr<FEMethod> null;
1068  CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), feAxiatorLhs,
1069  null, null);
1070  CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(), feAxiatorRhs,
1071  null, null);
1072  CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), feThermalLhs,
1073  null, null);
1074  CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(), feThermalRhs,
1075  null, null);
1076 
1077  CHKERR create_post_process_element();
1078  uXScatter = scatter_create(D, 0);
1079  uYScatter = scatter_create(D, 1);
1080  if (SPACE_DIM == 3)
1081  uZScatter = scatter_create(D, 2);
1082 
1083  auto solver = pipeline_mng->createTSIM();
1084 
1085  CHKERR TSSetSolution(solver, D);
1086  CHKERR set_section_monitor(solver);
1087  CHKERR set_time_monitor(dm, solver);
1088  CHKERR TSSetSolution(solver, D);
1089  CHKERR TSSetFromOptions(solver);
1090  CHKERR TSSetUp(solver);
1091  CHKERR TSSolve(solver, NULL);
1092 
1093  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
1094  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
1095  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
1096 
1098 }
1099 //! [Solve]
1100 
1101 static char help[] = "...\n\n";
1102 
1103 int main(int argc, char *argv[]) {
1104 
1105  // Initialisation of MoFEM/PETSc and MOAB data structures
1106  const char param_file[] = "param_file.petsc";
1107  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
1108 
1109  // Add logging channel for example
1110  auto core_log = logging::core::get();
1111  core_log->add_sink(
1113  LogManager::setLog("EXAMPLE");
1114  MOFEM_LOG_TAG("EXAMPLE", "example");
1115 
1116  try {
1117 
1118  //! [Register MoFEM discrete manager in PETSc]
1119  DMType dm_name = "DMMOFEM";
1120  CHKERR DMRegister_MoFEM(dm_name);
1121  //! [Register MoFEM discrete manager in PETSc
1122 
1123  //! [Create MoAB]
1124  moab::Core mb_instance; ///< mesh database
1125  moab::Interface &moab = mb_instance; ///< mesh database interface
1126  //! [Create MoAB]
1127 
1128  //! [Create MoFEM]
1129  MoFEM::Core core(moab); ///< finite element database
1130  MoFEM::Interface &m_field = core; ///< finite element database insterface
1131  //! [Create MoFEM]
1132 
1133  //! [Load mesh]
1134  Simple *simple = m_field.getInterface<Simple>();
1135  CHKERR simple->getOptions();
1136  CHKERR simple->loadFile();
1137  //! [Load mesh]
1138 
1139  //! [Example]
1140  Example ex(m_field);
1141  CHKERR ex.runProblem();
1142  //! [Example]
1143  }
1144  CATCH_ERRORS;
1145 
1147 }
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
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
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
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:79
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ H1
continuous field
Definition: definitions.h:98
@ HCURL
field with continuous tangents
Definition: definitions.h:99
@ HDIV
field with continuous normal traction
Definition: definitions.h:100
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
#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
constexpr auto t_kd
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1, 1 > OpBaseTimesScalarField
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:758
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:481
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:811
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:311
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:342
#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
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
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:994
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:1096
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
const double D
diffusivity
double A
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
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: plastic.cpp:74
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
Definition: plastic.cpp:95
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Body force]
Definition: plastic.cpp:72
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
Definition: plastic.cpp:97
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: plastic.cpp:93
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:41
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
PipelineManager::FaceEle DomainEle
BcTempFun(double v, FEMethod &fe)
double operator()(const double, const double, const double)
[Example]
std::vector< BcTempFun > bcTemperatureFunctions
boost::shared_ptr< PlasticThermalOps::CommonData > commonPlasticDataPtr
MoFEMErrorCode tsSolve()
boost::shared_ptr< DomainEle > feThermalRhs
boost::shared_ptr< DomainEle > feThermalLhs
MoFEMErrorCode createCommonData()
Example(MoFEM::Interface &m_field)
MoFEMErrorCode OPs()
MoFEMErrorCode runProblem()
MoFEMErrorCode setupProblem()
MoFEMErrorCode bC()
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
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
structure for User Loop Methods on finite elements
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 vector field for H-div approximation.
Calculate divergence of vector field.
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.
transform Hcurl base fluxes from reference element to physical triangle
Make Hdiv space from Hcurl space in 2d.
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]
double young_modulus
double heat_conductivity
int main(int argc, char *argv[])
double omega_inf
double amplitude_cycle
EntitiesFieldData::EntData EntData
double Qinf
static char help[]
[Solve]
double amplitude_shift
double phase_shift
double rho
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
#define EXECUTABLE_DIMENSION
constexpr int SPACE_DIM
double visH
double poisson_ratio
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< SPACE_DIM > OpHdivT
Integrate Lhs div of base of flux time base of temperature (FLUX x T) and transpose of it,...
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 3, 3, 1 > OpHdivFlux
Integrating Rhs flux base (1/k) flux (FLUX)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpBaseDotT
Integrate Rhs base of temerature time heat capacity times heat rate (T)
int number_of_cycles_in_total_time
double fraction_of_dissipation
double scale
long double hardening(long double tau, double temp)
double H
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpHeatSource
double omega_0
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpCapacity
Integrate Lhs base of temerature times (heat capacity) times base of temperature (T x T)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 1, 2 > OpHDivTemp
Integrate Rhs div flux base times temperature (T)
double heat_capacity
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, 3 > OpHdivHdiv
Integrate Lhs base of flux (1/k) base of flux (FLUX x FLUX)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Body force]
long double hardening_dtau(long double tau, double temp)
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
int order
double b_iso
double sigmaY
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpNormalMixVecTimesScalar< SPACE_DIM > OpTemperatureBC
long double hardening_dtemp(long double tau, double temp)
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
double cn
OpBaseDotT OpBaseDivFlux
Integrate Rhs base of temerature times divergenc of flux (T)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, SPACE_DIM > OpBodyForce
[Body force]
constexpr EntityType boundary_ent
double omega_h