v0.14.0
thermo_elastic.cpp
Go to the documentation of this file.
1 /**
2  * \file thermo_elastic.cpp
3  * \example thermo_elastic.cpp
4  *
5  * Thermo plasticity
6  *
7  */
8 
9 #ifndef EXECUTABLE_DIMENSION
10 #define EXECUTABLE_DIMENSION 3
11 #endif
12 
13 #include <MoFEM.hpp>
14 
15 using namespace MoFEM;
16 
17 constexpr int SPACE_DIM =
18  EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
19 
20 using EntData = EntitiesFieldData::EntData;
21 using DomainEle = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::DomainEle;
22 using DomainEleOp = DomainEle::UserDataOperator;
23 using BoundaryEle =
24  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
25 using BoundaryEleOp = BoundaryEle::UserDataOperator;
26 using PostProcEle = PostProcBrokenMeshInMoab<DomainEle>;
27 
28 using AssemblyDomainEleOp =
29  FormsIntegrators<DomainEleOp>::Assembly<PETSC>::OpBase;
30 
31 //! [Linear elastic problem]
32 using OpKCauchy = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
33  GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM,
34  0>; //< Elastic stiffness matrix
35 using OpInternalForceCauchy =
36  FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
37  GAUSS>::OpGradTimesSymTensor<1, SPACE_DIM,
38  SPACE_DIM>; //< Elastic internal forces
39 //! [Linear elastic problem]
40 
41 //! [Thermal problem]
42 /**
43  * @brief Integrate Lhs base of flux (1/k) base of flux (FLUX x FLUX)
44  *
45  */
46 using OpHdivHdiv = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
47  GAUSS>::OpMass<3, SPACE_DIM>;
48 
49 /**
50  * @brief Integrate Lhs div of base of flux time base of temperature (FLUX x T)
51  * and transpose of it, i.e. (T x FLAX)
52  *
53  */
54 using OpHdivT = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
55  GAUSS>::OpMixDivTimesScalar<SPACE_DIM>;
56 
57 /**
58  * @brief Integrate Lhs base of temperature times (heat capacity) times base of
59  * temperature (T x T)
60  *
61  */
62 using OpCapacity = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
63  GAUSS>::OpMass<1, 1>;
64 
65 /**
66  * @brief Integrating Rhs flux base (1/k) flux (FLUX)
67  */
68 using OpHdivFlux = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
69  GAUSS>::OpBaseTimesVector<3, SPACE_DIM, 1>;
70 
71 /**
72  * @brief Integrate Rhs div flux base times temperature (T)
73  *
74  */
75 using OpHDivTemp = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
76  GAUSS>::OpMixDivTimesU<3, 1, SPACE_DIM>;
77 
78 /**
79  * @brief Integrate Rhs base of temperature time heat capacity times heat rate
80  * (T)
81  *
82  */
83 using OpBaseDotT = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
84  GAUSS>::OpBaseTimesScalar<1>;
85 
86 /**
87  * @brief Integrate Rhs base of temperature times divergence of flux (T)
88  *
89  */
90 using OpBaseDivFlux = OpBaseDotT;
91 
92 //! [Thermal problem]
93 
94 //! [Body and heat source]
95 using DomainNaturalBCRhs =
96  NaturalBC<DomainEleOp>::Assembly<PETSC>::LinearForm<GAUSS>;
97 using OpBodyForce =
98  DomainNaturalBCRhs::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, SPACE_DIM>;
99 using OpHeatSource =
100  DomainNaturalBCRhs::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 1>;
101 using DomainNaturalBCLhs =
102  NaturalBC<DomainEleOp>::Assembly<PETSC>::BiLinearForm<GAUSS>;
103 //! [Body and heat source]
104 
105 //! [Natural boundary conditions]
106 using BoundaryNaturalBC =
107  NaturalBC<BoundaryEleOp>::Assembly<PETSC>::LinearForm<GAUSS>;
108 using OpForce = BoundaryNaturalBC::OpFlux<NaturalForceMeshsets, 1, SPACE_DIM>;
109 using OpTemperatureBC =
110  BoundaryNaturalBC::OpFlux<NaturalTemperatureMeshsets, 3, SPACE_DIM>;
111 //! [Natural boundary conditions]
112 
113 //! [Essential boundary conditions (Least square approach)]
114 using OpEssentialFluxRhs =
115  EssentialBC<BoundaryEleOp>::Assembly<PETSC>::LinearForm<
116  GAUSS>::OpEssentialRhs<HeatFluxCubitBcData, 3, SPACE_DIM>;
117 using OpEssentialFluxLhs =
118  EssentialBC<BoundaryEleOp>::Assembly<PETSC>::BiLinearForm<
119  GAUSS>::OpEssentialLhs<HeatFluxCubitBcData, 3, SPACE_DIM>;
120 //! [Essential boundary conditions (Least square approach)]
121 
122 double default_young_modulus = 1;
123 double default_poisson_ratio = 0.25;
124 double ref_temp = 0.0;
125 double init_temp = 0.0;
126 
127 double default_coeff_expansion = 1;
128 double default_heat_conductivity =
129  1; // Force / (time temperature ) or Power /
130  // (length temperature) // Time unit is hour. force unit kN
131 double default_heat_capacity = 1; // length^2/(time^2 temperature) // length is
132  // millimeter time is hour
133 
134 int order = 2; //< default approximation order
135 
136 #include <ThermoElasticOps.hpp> //< additional coupling opearyors
137 using namespace ThermoElasticOps; //< name space of coupling operators
138 
139 using OpSetTemperatureRhs =
140  DomainNaturalBCRhs::OpFlux<SetTargetTemperature, 1, 1>;
141 using OpSetTemperatureLhs =
142  DomainNaturalBCLhs::OpFlux<SetTargetTemperature, 1, 1>;
143 
144 auto save_range = [](moab::Interface &moab, const std::string name,
145  const Range r) {
146  MoFEMFunctionBegin;
147  auto out_meshset = get_temp_meshset_ptr(moab);
148  CHKERR moab.add_entities(*out_meshset, r);
149  CHKERR moab.write_file(name.c_str(), "VTK", "", out_meshset->get_ptr(), 1);
150  MoFEMFunctionReturn(0);
151 };
152 
153 struct ThermoElasticProblem {
154 
155  ThermoElasticProblem(MoFEM::Interface &m_field) : mField(m_field) {}
156 
157  MoFEMErrorCode runProblem();
158 
159 private:
160  MoFEM::Interface &mField;
161 
162  PetscBool doEvalField;
163  std::array<double, SPACE_DIM> fieldEvalCoords;
164  boost::shared_ptr<FieldEvaluatorInterface::SetPtsData> fieldEvalData;
165 
166  boost::shared_ptr<VectorDouble> scalarFieldPtr;
167  boost::shared_ptr<MatrixDouble> vectorFieldPtr;
168  boost::shared_ptr<MatrixDouble> tensorFieldPtr;
169 
170  MoFEMErrorCode setupProblem(); ///< add fields
171  MoFEMErrorCode createCommonData(); //< read global data from command line
172  MoFEMErrorCode bC(); //< read boundary conditions
173  MoFEMErrorCode OPs(); //< add operators to pipeline
174  MoFEMErrorCode tsSolve(); //< time solver
175 
176  struct BlockedParameters
177  : public boost::enable_shared_from_this<BlockedParameters> {
178  MatrixDouble mD;
179  double coeffExpansion;
180  double heatConductivity;
181  double heatCapacity;
182 
183  inline auto getDPtr() {
184  return boost::shared_ptr<MatrixDouble>(shared_from_this(), &mD);
185  }
186 
187  inline auto getCoeffExpansionPtr() {
188  return boost::shared_ptr<double>(shared_from_this(), &coeffExpansion);
189  }
190 
191  inline auto getHeatConductivityPtr() {
192  return boost::shared_ptr<double>(shared_from_this(), &heatConductivity);
193  }
194 
195  inline auto getHeatCapacityPtr() {
196  return boost::shared_ptr<double>(shared_from_this(), &heatCapacity);
197  }
198  };
199 
200  MoFEMErrorCode addMatBlockOps(
201  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
202  std::string block_elastic_name, std::string block_thermal_name,
203  boost::shared_ptr<BlockedParameters> blockedParamsPtr, Sev sev);
204 };
205 
206 MoFEMErrorCode ThermoElasticProblem::addMatBlockOps(
207  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
208  std::string block_elastic_name, std::string block_thermal_name,
209  boost::shared_ptr<BlockedParameters> blockedParamsPtr, Sev sev) {
210  MoFEMFunctionBegin;
211 
212  struct OpMatElasticBlocks : public DomainEleOp {
213  OpMatElasticBlocks(boost::shared_ptr<MatrixDouble> m, double bulk_modulus_K,
214  double shear_modulus_G, MoFEM::Interface &m_field,
215  Sev sev,
216  std::vector<const CubitMeshSets *> meshset_vec_ptr)
217  : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), matDPtr(m),
218  bulkModulusKDefault(bulk_modulus_K),
219  shearModulusGDefault(shear_modulus_G) {
220  CHK_THROW_MESSAGE(extractElasticBlockData(m_field, meshset_vec_ptr, sev),
221  "Can not get data from block");
222  }
223 
224  MoFEMErrorCode doWork(int side, EntityType type,
225  EntitiesFieldData::EntData &data) {
226  MoFEMFunctionBegin;
227 
228  for (auto &b : blockData) {
229 
230  if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
231  CHKERR getMatDPtr(matDPtr, b.bulkModulusK, b.shearModulusG);
232  MoFEMFunctionReturnHot(0);
233  }
234  }
235 
236  CHKERR getMatDPtr(matDPtr, bulkModulusKDefault, shearModulusGDefault);
237  MoFEMFunctionReturn(0);
238  }
239 
240  private:
241  boost::shared_ptr<MatrixDouble> matDPtr;
242 
243  struct BlockData {
244  double bulkModulusK;
245  double shearModulusG;
246  Range blockEnts;
247  };
248 
249  double bulkModulusKDefault;
250  double shearModulusGDefault;
251  std::vector<BlockData> blockData;
252 
253  MoFEMErrorCode
254  extractElasticBlockData(MoFEM::Interface &m_field,
255  std::vector<const CubitMeshSets *> meshset_vec_ptr,
256  Sev sev) {
257  MoFEMFunctionBegin;
258 
259  for (auto m : meshset_vec_ptr) {
260  MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Elastic Block") << *m;
261  std::vector<double> block_data;
262  CHKERR m->getAttributes(block_data);
263  if (block_data.size() < 2) {
264  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
265  "Expected that block has two attributes");
266  }
267  auto get_block_ents = [&]() {
268  Range ents;
269  CHKERR
270  m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
271  return ents;
272  };
273 
274  double young_modulus = block_data[0];
275  double poisson_ratio = block_data[1];
276  double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
277  double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
278 
279  MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Elastic Block")
280  << m->getName() << ": E = " << young_modulus
281  << " nu = " << poisson_ratio;
282 
283  blockData.push_back(
284  {bulk_modulus_K, shear_modulus_G, get_block_ents()});
285  }
286  MOFEM_LOG_CHANNEL("WORLD");
287  MoFEMFunctionReturn(0);
288  }
289 
290  MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
291  double bulk_modulus_K, double shear_modulus_G) {
292  MoFEMFunctionBegin;
293  //! [Calculate elasticity tensor]
294  auto set_material_stiffness = [&]() {
295  FTensor::Index<'i', SPACE_DIM> i;
296  FTensor::Index<'j', SPACE_DIM> j;
297  FTensor::Index<'k', SPACE_DIM> k;
298  FTensor::Index<'l', SPACE_DIM> l;
299  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
300  double A = (SPACE_DIM == 2)
301  ? 2 * shear_modulus_G /
302  (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
303  : 1;
304  auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
305  t_D(i, j, k, l) =
306  2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
307  A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) * t_kd(i, j) *
308  t_kd(k, l);
309  };
310  //! [Calculate elasticity tensor]
311  constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
312  mat_D_ptr->resize(size_symm * size_symm, 1);
313  set_material_stiffness();
314  MoFEMFunctionReturn(0);
315  }
316  };
317 
318  double default_bulk_modulus_K =
319  default_young_modulus / (3 * (1 - 2 * default_poisson_ratio));
320  double default_shear_modulus_G =
321  default_young_modulus / (2 * (1 + default_poisson_ratio));
322 
323  pipeline.push_back(new OpMatElasticBlocks(
324  blockedParamsPtr->getDPtr(), default_bulk_modulus_K,
325  default_bulk_modulus_K, mField, sev,
326 
327  // Get blockset using regular expression
328  mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
329 
330  (boost::format("%s(.*)") % block_elastic_name).str()
331 
332  ))
333 
334  ));
335 
336  struct OpMatThermalBlocks : public DomainEleOp {
337  OpMatThermalBlocks(boost::shared_ptr<double> expansion_ptr,
338  boost::shared_ptr<double> conductivity_ptr,
339  boost::shared_ptr<double> capacity_ptr,
340  MoFEM::Interface &m_field, Sev sev,
341  std::vector<const CubitMeshSets *> meshset_vec_ptr)
342  : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE),
343  expansionPtr(expansion_ptr), conductivityPtr(conductivity_ptr),
344  capacityPtr(capacity_ptr) {
345  CHK_THROW_MESSAGE(extractThermallockData(m_field, meshset_vec_ptr, sev),
346  "Can not get data from block");
347  }
348 
349  MoFEMErrorCode doWork(int side, EntityType type,
350  EntitiesFieldData::EntData &data) {
351  MoFEMFunctionBegin;
352 
353  for (auto &b : blockData) {
354 
355  if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
356  *expansionPtr = b.expansion;
357  *conductivityPtr = b.conductivity;
358  *capacityPtr = b.capacity;
359  MoFEMFunctionReturnHot(0);
360  }
361  }
362 
363  *expansionPtr = default_coeff_expansion;
364  *conductivityPtr = default_heat_conductivity;
365  *capacityPtr = default_heat_capacity;
366 
367  MoFEMFunctionReturn(0);
368  }
369 
370  private:
371  struct BlockData {
372  double expansion;
373  double conductivity;
374  double capacity;
375  Range blockEnts;
376  };
377 
378  std::vector<BlockData> blockData;
379 
380  MoFEMErrorCode
381  extractThermallockData(MoFEM::Interface &m_field,
382  std::vector<const CubitMeshSets *> meshset_vec_ptr,
383  Sev sev) {
384  MoFEMFunctionBegin;
385 
386  for (auto m : meshset_vec_ptr) {
387  MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Thermal Block") << *m;
388  std::vector<double> block_data;
389  CHKERR m->getAttributes(block_data);
390  if (block_data.size() < 3) {
391  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
392  "Expected that block has two attributes");
393  }
394  auto get_block_ents = [&]() {
395  Range ents;
396  CHKERR
397  m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
398  return ents;
399  };
400 
401  MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Thermal Block")
402  << m->getName() << ": expansion = " << block_data[0]
403  << " conductivity = " << block_data[1] << " capacity "
404  << block_data[2];
405 
406  blockData.push_back(
407  {block_data[0], block_data[1], block_data[2], get_block_ents()});
408  }
409  MOFEM_LOG_CHANNEL("WORLD");
410  MoFEMFunctionReturn(0);
411  }
412 
413  boost::shared_ptr<double> expansionPtr;
414  boost::shared_ptr<double> conductivityPtr;
415  boost::shared_ptr<double> capacityPtr;
416  };
417 
418  pipeline.push_back(new OpMatThermalBlocks(
419  blockedParamsPtr->getCoeffExpansionPtr(),
420  blockedParamsPtr->getHeatConductivityPtr(),
421  blockedParamsPtr->getHeatCapacityPtr(), mField, sev,
422 
423  // Get blockset using regular expression
424  mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
425 
426  (boost::format("%s(.*)") % block_thermal_name).str()
427 
428  ))
429 
430  ));
431 
432  MoFEMFunctionReturn(0);
433 }
434 
435 //! [Run problem]
436 MoFEMErrorCode ThermoElasticProblem::runProblem() {
437  MoFEMFunctionBegin;
438  CHKERR setupProblem();
439  CHKERR createCommonData();
440  CHKERR bC();
441  CHKERR OPs();
442  CHKERR tsSolve();
443  MoFEMFunctionReturn(0);
444 }
445 //! [Run problem]
446 
447 //! [Set up problem]
448 MoFEMErrorCode ThermoElasticProblem::setupProblem() {
449  MoFEMFunctionBegin;
450  Simple *simple = mField.getInterface<Simple>();
451  // Add field
452  constexpr FieldApproximationBase base = AINSWORTH_LEGENDRE_BASE;
453  // Mechanical fields
454  CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
455  CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
456  // Temperature
457  constexpr auto flux_space = (SPACE_DIM == 2) ? HCURL : HDIV;
458  CHKERR simple->addDomainField("T", L2, AINSWORTH_LEGENDRE_BASE, 1);
459  CHKERR simple->addDomainField("FLUX", flux_space, DEMKOWICZ_JACOBI_BASE, 1);
460  CHKERR simple->addBoundaryField("FLUX", flux_space, DEMKOWICZ_JACOBI_BASE, 1);
461 
462  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
463  CHKERR simple->setFieldOrder("U", order + 1);
464  CHKERR simple->setFieldOrder("FLUX", order + 1);
465  CHKERR simple->setFieldOrder("T", order);
466  CHKERR simple->setUp();
467 
468  int coords_dim = SPACE_DIM;
469  CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
470  fieldEvalCoords.data(), &coords_dim,
471  &doEvalField);
472 
473  scalarFieldPtr = boost::make_shared<VectorDouble>();
474  vectorFieldPtr = boost::make_shared<MatrixDouble>();
475  tensorFieldPtr = boost::make_shared<MatrixDouble>();
476 
477  if (doEvalField) {
478  fieldEvalData =
479  mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
480 
481  if constexpr (SPACE_DIM == 3) {
482  CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree3D(
483  fieldEvalData, simple->getDomainFEName());
484  } else {
485  CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree2D(
486  fieldEvalData, simple->getDomainFEName());
487  }
488 
489  fieldEvalData->setEvalPoints(fieldEvalCoords.data(), 1);
490  auto no_rule = [](int, int, int) { return -1; };
491 
492  auto field_eval_fe_ptr = fieldEvalData->feMethodPtr.lock();
493  field_eval_fe_ptr->getRuleHook = no_rule;
494 
495  field_eval_fe_ptr->getOpPtrVector().push_back(
496  new OpCalculateScalarFieldValues("T", scalarFieldPtr));
497  field_eval_fe_ptr->getOpPtrVector().push_back(
498  new OpCalculateVectorFieldValues<SPACE_DIM>("U", vectorFieldPtr));
499  field_eval_fe_ptr->getOpPtrVector().push_back(
500  new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
501  "U", tensorFieldPtr));
502  }
503 
504  MoFEMFunctionReturn(0);
505 }
506 //! [Set up problem]
507 
508 //! [Create common data]
509 MoFEMErrorCode ThermoElasticProblem::createCommonData() {
510  MoFEMFunctionBegin;
511 
512  auto get_command_line_parameters = [&]() {
513  MoFEMFunctionBegin;
514  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
515  &default_young_modulus, PETSC_NULL);
516  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
517  &default_poisson_ratio, PETSC_NULL);
518  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-coeff_expansion",
519  &default_coeff_expansion, PETSC_NULL);
520  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-ref_temp", &ref_temp,
521  PETSC_NULL);
522  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-init_temp", &init_temp,
523  PETSC_NULL);
524 
525  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-capacity",
526  &default_heat_capacity, PETSC_NULL);
527  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-conductivity",
528  &default_heat_conductivity, PETSC_NULL);
529 
530  MOFEM_LOG("ThermoElastic", Sev::inform)
531  << "Young modulus " << default_young_modulus;
532  MOFEM_LOG("ThermoElastic", Sev::inform)
533  << "Poisson ratio " << default_poisson_ratio;
534  MOFEM_LOG("ThermoElastic", Sev::inform)
535  << "Coeff of expansion " << default_coeff_expansion;
536  MOFEM_LOG("ThermoElastic", Sev::inform)
537  << "Capacity " << default_heat_capacity;
538  MOFEM_LOG("ThermoElastic", Sev::inform)
539  << "Heat conductivity " << default_heat_conductivity;
540 
541  MOFEM_LOG("ThermoElastic", Sev::inform)
542  << "Reference temperature " << ref_temp;
543  MOFEM_LOG("ThermoElastic", Sev::inform)
544  << "Initial temperature " << init_temp;
545 
546  MoFEMFunctionReturn(0);
547  };
548 
549  CHKERR get_command_line_parameters();
550 
551  MoFEMFunctionReturn(0);
552 }
553 //! [Create common data]
554 
555 //! [Boundary condition]
556 MoFEMErrorCode ThermoElasticProblem::bC() {
557  MoFEMFunctionBegin;
558 
559  MOFEM_LOG("SYNC", Sev::noisy) << "bC";
560  MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::noisy);
561 
562  auto simple = mField.getInterface<Simple>();
563  auto bc_mng = mField.getInterface<BcManager>();
564 
565 
566 
567  auto get_skin = [&]() {
568  Range body_ents;
569  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
570  Skinner skin(&mField.get_moab());
571  Range skin_ents;
572  CHKERR skin.find_skin(0, body_ents, false, skin_ents);
573  return skin_ents;
574  };
575 
576  auto filter_flux_blocks = [&](auto skin) {
577  auto remove_cubit_blocks = [&](auto c) {
578  MoFEMFunctionBegin;
579  for (auto m :
580 
581  mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(c)
582 
583  ) {
584  Range ents;
585  CHKERR mField.get_moab().get_entities_by_dimension(
586  m->getMeshset(), SPACE_DIM - 1, ents, true);
587  skin = subtract(skin, ents);
588  }
589  MoFEMFunctionReturn(0);
590  };
591 
592  auto remove_named_blocks = [&](auto n) {
593  MoFEMFunctionBegin;
594  for (auto m : mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
595  std::regex(
596 
597  (boost::format("%s(.*)") % n).str()
598 
599  ))
600 
601  ) {
602  Range ents;
603  CHKERR mField.get_moab().get_entities_by_dimension(
604  m->getMeshset(), SPACE_DIM - 1, ents, true);
605  skin = subtract(skin, ents);
606  }
607  MoFEMFunctionReturn(0);
608  };
609 
610  CHK_THROW_MESSAGE(remove_cubit_blocks(NODESET | TEMPERATURESET),
611  "remove_cubit_blocks");
612  CHK_THROW_MESSAGE(remove_cubit_blocks(SIDESET | HEATFLUXSET),
613  "remove_cubit_blocks");
614  CHK_THROW_MESSAGE(remove_named_blocks("TEMPERATURE"),
615  "remove_named_blocks");
616  CHK_THROW_MESSAGE(remove_named_blocks("HEATFLUX"), "remove_named_blocks");
617 
618  return skin;
619  };
620 
621  auto filter_true_skin = [&](auto skin) {
622  Range boundary_ents;
623  ParallelComm *pcomm =
624  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
625  CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
626  PSTATUS_NOT, -1, &boundary_ents);
627  return boundary_ents;
628  };
629 
630  auto remove_flux_ents = filter_true_skin(filter_flux_blocks(get_skin()));
631 
632  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
633  remove_flux_ents);
634 
635  MOFEM_LOG("SYNC", Sev::noisy) << remove_flux_ents << endl;
636  MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::noisy);
637 
638 #ifndef NDEBUG
639 
640  CHKERR save_range(
641  mField.get_moab(),
642  (boost::format("flux_remove_%d.vtk") % mField.get_comm_rank()).str(),
643  remove_flux_ents);
644 
645 #endif
646 
647  CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
648  simple->getProblemName(), "FLUX", remove_flux_ents);
649 
650  auto set_init_temp = [](boost::shared_ptr<FieldEntity> field_entity_ptr) {
651  field_entity_ptr->getEntFieldData()[0] = init_temp;
652  return 0;
653  };
654  CHKERR mField.getInterface<FieldBlas>()->fieldLambdaOnEntities(set_init_temp,
655  "T");
656 
657  CHKERR bc_mng->removeBlockDOFsOnEntities<DisplacementCubitBcData>(
658  simple->getProblemName(), "U");
659  CHKERR bc_mng->pushMarkDOFsOnEntities<HeatFluxCubitBcData>(
660  simple->getProblemName(), "FLUX", false);
661 
662  MoFEMFunctionReturn(0);
663 }
664 //! [Boundary condition]
665 
666 //! [Push operators to pipeline]
667 MoFEMErrorCode ThermoElasticProblem::OPs() {
668  MoFEMFunctionBegin;
669 
670  MOFEM_LOG("SYNC", Sev::noisy) << "OPs";
671  MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::noisy);
672 
673  auto pipeline_mng = mField.getInterface<PipelineManager>();
674  auto simple = mField.getInterface<Simple>();
675  auto bc_mng = mField.getInterface<BcManager>();
676 
677  auto boundary_marker =
678  bc_mng->getMergedBlocksMarker(vector<string>{"HEATFLUX"});
679 
680  auto integration_rule = [](int, int, int approx_order) {
681  return 2 * approx_order;
682  };
683  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
684  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
685  CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
686  CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
687 
688  auto block_params = boost::make_shared<BlockedParameters>();
689  auto mDPtr = block_params->getDPtr();
690  auto coeff_expansion_ptr = block_params->getCoeffExpansionPtr();
691  auto heat_conductivity_ptr = block_params->getHeatConductivityPtr();
692  auto heat_capacity_ptr = block_params->getHeatCapacityPtr();
693 
694  // Default time scaling of BCs and sources, from command line arguments
695  auto time_scale = boost::make_shared<TimeScale>();
696 
697  // Files which define scaling for separate variables, if provided
698  auto time_bodyforce_scaling =
699  boost::make_shared<TimeScale>("bodyforce_scale.txt");
700  auto time_heatsource_scaling =
701  boost::make_shared<TimeScale>("heatsource_scale.txt");
702  auto time_temperature_scaling =
703  boost::make_shared<TimeScale>("temperature_bc_scale.txt");
704  auto time_displacement_scaling =
705  boost::make_shared<TimeScale>("displacement_bc_scale.txt");
706  auto time_flux_scaling = boost::make_shared<TimeScale>("flux_bc_scale.txt");
707  auto time_force_scaling =
708  boost::make_shared<TimeScale>("force_bc_scale.txt");
709 
710  auto add_domain_rhs_ops = [&](auto &pipeline) {
711  MoFEMFunctionBegin;
712  CHKERR addMatBlockOps(pipeline, "MAT_ELASTIC", "MAT_THERMAL", block_params,
713  Sev::inform);
714  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pipeline, {H1, HDIV});
715 
716  auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
717  auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
718  auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
719 
720  auto vec_temp_ptr = boost::make_shared<VectorDouble>();
721  auto vec_temp_dot_ptr = boost::make_shared<VectorDouble>();
722  auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
723  auto vec_temp_div_ptr = boost::make_shared<VectorDouble>();
724 
725  pipeline.push_back(new OpCalculateScalarFieldValues("T", vec_temp_ptr));
726  pipeline.push_back(
727  new OpCalculateScalarFieldValuesDot("T", vec_temp_dot_ptr));
728  pipeline.push_back(new OpCalculateHdivVectorDivergence<3, SPACE_DIM>(
729  "FLUX", vec_temp_div_ptr));
730  pipeline.push_back(
731  new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
732 
733  pipeline.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
734  "U", mat_grad_ptr));
735  pipeline.push_back(
736  new OpSymmetrizeTensor<SPACE_DIM>(mat_grad_ptr, mat_strain_ptr));
737  pipeline.push_back(new OpStressThermal(mat_strain_ptr, vec_temp_ptr, mDPtr,
738  coeff_expansion_ptr,
739  mat_stress_ptr));
740 
741  pipeline.push_back(new OpSetBc("FLUX", true, boundary_marker));
742 
743  pipeline.push_back(new OpInternalForceCauchy(
744  "U", mat_stress_ptr,
745  [](double, double, double) constexpr { return 1; }));
746 
747  auto resistance = [heat_conductivity_ptr](const double, const double,
748  const double) {
749  return (1. / (*heat_conductivity_ptr));
750  };
751  // negative value is consequence of symmetric system
752  auto capacity = [heat_capacity_ptr](const double, const double,
753  const double) {
754  return -(*heat_capacity_ptr);
755  };
756  auto unity = [](const double, const double, const double) constexpr {
757  return -1.;
758  };
759  pipeline.push_back(new OpHdivFlux("FLUX", mat_flux_ptr, resistance));
760  pipeline.push_back(new OpHDivTemp("FLUX", vec_temp_ptr, unity));
761  pipeline.push_back(new OpBaseDivFlux("T", vec_temp_div_ptr, unity));
762  pipeline.push_back(new OpBaseDotT("T", vec_temp_dot_ptr, capacity));
763 
764  pipeline.push_back(new OpUnSetBc("FLUX"));
765 
766  // Set source terms
767  CHKERR DomainNaturalBCRhs::AddFluxToPipeline<OpHeatSource>::add(
768  pipeline, mField, "T", {time_scale, time_heatsource_scaling},
769  "HEAT_SOURCE", Sev::inform);
770  CHKERR DomainNaturalBCRhs::AddFluxToPipeline<OpBodyForce>::add(
771  pipeline, mField, "U", {time_scale, time_bodyforce_scaling},
772  "BODY_FORCE", Sev::inform);
773  CHKERR DomainNaturalBCRhs::AddFluxToPipeline<OpSetTemperatureRhs>::add(
774  pipeline, mField, "T", vec_temp_ptr, "SETTEMP", Sev::inform);
775 
776  MoFEMFunctionReturn(0);
777  };
778 
779  auto add_domain_lhs_ops = [&](auto &pipeline) {
780  MoFEMFunctionBegin;
781  CHKERR addMatBlockOps(pipeline, "MAT_ELASTIC", "MAT_THERMAL", block_params,
782  Sev::verbose);
783  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pipeline, {H1, HDIV});
784 
785  pipeline.push_back(new OpSetBc("FLUX", true, boundary_marker));
786 
787  pipeline.push_back(new OpKCauchy("U", "U", mDPtr));
788  pipeline.push_back(new ThermoElasticOps::OpKCauchyThermoElasticity(
789  "U", "T", mDPtr, coeff_expansion_ptr));
790 
791  auto resistance = [heat_conductivity_ptr](const double, const double,
792  const double) {
793  return (1. / (*heat_conductivity_ptr));
794  };
795  auto capacity = [heat_capacity_ptr](const double, const double,
796  const double) {
797  return -(*heat_capacity_ptr);
798  };
799  pipeline.push_back(new OpHdivHdiv("FLUX", "FLUX", resistance));
800  pipeline.push_back(new OpHdivT(
801  "FLUX", "T", []() constexpr { return -1; }, true));
802 
803  auto op_capacity = new OpCapacity("T", "T", capacity);
804  op_capacity->feScalingFun = [](const FEMethod *fe_ptr) {
805  return fe_ptr->ts_a;
806  };
807  pipeline.push_back(op_capacity);
808 
809  pipeline.push_back(new OpUnSetBc("FLUX"));
810 
811  auto vec_temp_ptr = boost::make_shared<VectorDouble>();
812  pipeline.push_back(new OpCalculateScalarFieldValues("T", vec_temp_ptr));
813  CHKERR DomainNaturalBCLhs::AddFluxToPipeline<OpSetTemperatureLhs>::add(
814  pipeline, mField, "T", vec_temp_ptr, "SETTEMP", Sev::verbose);
815 
816  MoFEMFunctionReturn(0);
817  };
818 
819  auto add_boundary_rhs_ops = [&](auto &pipeline) {
820  MoFEMFunctionBegin;
821 
822  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pipeline, {HDIV});
823 
824  pipeline.push_back(new OpSetBc("FLUX", true, boundary_marker));
825 
826  // Set BCs using the least squares imposition
827  CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpForce>::add(
828  pipeline_mng->getOpBoundaryRhsPipeline(), mField, "U",
829  {time_scale, time_force_scaling}, "FORCE", Sev::inform);
830 
831  CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpTemperatureBC>::add(
832  pipeline_mng->getOpBoundaryRhsPipeline(), mField, "FLUX",
833  {time_scale, time_temperature_scaling}, "TEMPERATURE", Sev::inform);
834 
835  pipeline.push_back(new OpUnSetBc("FLUX"));
836 
837  auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
838  pipeline.push_back(
839  new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
840  CHKERR EssentialBC<BoundaryEleOp>::Assembly<PETSC>::LinearForm<GAUSS>::
841  AddEssentialToPipeline<OpEssentialFluxRhs>::add(
842  mField, pipeline, simple->getProblemName(), "FLUX", mat_flux_ptr,
843  {time_scale, time_flux_scaling});
844 
845  MoFEMFunctionReturn(0);
846  };
847 
848  auto add_boundary_lhs_ops = [&](auto &pipeline) {
849  MoFEMFunctionBegin;
850 
851  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pipeline, {HDIV});
852 
853  CHKERR EssentialBC<BoundaryEleOp>::Assembly<PETSC>::BiLinearForm<GAUSS>::
854  AddEssentialToPipeline<OpEssentialFluxLhs>::add(
855  mField, pipeline, simple->getProblemName(), "FLUX");
856 
857  MoFEMFunctionReturn(0);
858  };
859 
860  // Set BCs by eliminating degrees of freedom
861  auto get_bc_hook_rhs = [&]() {
862  EssentialPreProc<DisplacementCubitBcData> hook(
863  mField, pipeline_mng->getDomainRhsFE(),
864  {time_scale, time_displacement_scaling});
865  return hook;
866  };
867  auto get_bc_hook_lhs = [&]() {
868  EssentialPreProc<DisplacementCubitBcData> hook(
869  mField, pipeline_mng->getDomainLhsFE(),
870  {time_scale, time_displacement_scaling});
871  return hook;
872  };
873 
874  pipeline_mng->getDomainRhsFE()->preProcessHook = get_bc_hook_rhs();
875  pipeline_mng->getDomainLhsFE()->preProcessHook = get_bc_hook_lhs();
876 
877  CHKERR add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
878  CHKERR add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
879  CHKERR add_boundary_rhs_ops(pipeline_mng->getOpBoundaryRhsPipeline());
880  CHKERR add_boundary_lhs_ops(pipeline_mng->getOpBoundaryLhsPipeline());
881 
882  MoFEMFunctionReturn(0);
883 }
884 //! [Push operators to pipeline]
885 
886 //! [Solve]
887 MoFEMErrorCode ThermoElasticProblem::tsSolve() {
888  MoFEMFunctionBegin;
889 
890  Simple *simple = mField.getInterface<Simple>();
891  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
892  ISManager *is_manager = mField.getInterface<ISManager>();
893 
894  auto dm = simple->getDM();
895  auto solver = pipeline_mng->createTSIM();
896  auto snes_ctx_ptr = getDMSnesCtx(dm);
897 
898  auto set_section_monitor = [&](auto solver) {
899  MoFEMFunctionBegin;
900  SNES snes;
901  CHKERR TSGetSNES(solver, &snes);
902  CHKERR SNESMonitorSet(snes,
903  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
904  void *))MoFEMSNESMonitorFields,
905  (void *)(snes_ctx_ptr.get()), nullptr);
906  MoFEMFunctionReturn(0);
907  };
908 
909  auto create_post_process_element = [&]() {
910  auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
911 
912  auto block_params = boost::make_shared<BlockedParameters>();
913  auto mDPtr = block_params->getDPtr();
914  auto coeff_expansion_ptr = block_params->getCoeffExpansionPtr();
915 
916  CHKERR addMatBlockOps(post_proc_fe->getOpPtrVector(), "MAT_ELASTIC",
917  "MAT_THERMAL", block_params, Sev::verbose);
918 
919  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
920  post_proc_fe->getOpPtrVector(), {H1, HDIV});
921 
922  auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
923  auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
924  auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
925 
926  auto vec_temp_ptr = boost::make_shared<VectorDouble>();
927  auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
928 
929  post_proc_fe->getOpPtrVector().push_back(
930  new OpCalculateScalarFieldValues("T", vec_temp_ptr));
931  post_proc_fe->getOpPtrVector().push_back(
932  new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
933 
934  auto u_ptr = boost::make_shared<MatrixDouble>();
935  post_proc_fe->getOpPtrVector().push_back(
936  new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
937  post_proc_fe->getOpPtrVector().push_back(
938  new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>("U",
939  mat_grad_ptr));
940  post_proc_fe->getOpPtrVector().push_back(
941  new OpSymmetrizeTensor<SPACE_DIM>(mat_grad_ptr, mat_strain_ptr));
942  post_proc_fe->getOpPtrVector().push_back(
943  new OpStressThermal(mat_strain_ptr, vec_temp_ptr, mDPtr,
944  coeff_expansion_ptr, mat_stress_ptr));
945 
946  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
947 
948  post_proc_fe->getOpPtrVector().push_back(
949 
950  new OpPPMap(
951 
952  post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
953 
954  {{"T", vec_temp_ptr}},
955 
956  {{"U", u_ptr}, {"FLUX", mat_flux_ptr}},
957 
958  {},
959 
960  {{"STRAIN", mat_strain_ptr}, {"STRESS", mat_stress_ptr}}
961 
962  )
963 
964  );
965 
966  return post_proc_fe;
967  };
968 
969  auto monitor_ptr = boost::make_shared<FEMethod>();
970  auto post_proc_fe = create_post_process_element();
971 
972  auto set_time_monitor = [&](auto dm, auto solver) {
973  MoFEMFunctionBegin;
974  monitor_ptr->preProcessHook = [&]() {
975  MoFEMFunctionBegin;
976 
977  CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
978  post_proc_fe,
979  monitor_ptr->getCacheWeakPtr());
980  CHKERR post_proc_fe->writeFile(
981  "out_" + boost::lexical_cast<std::string>(monitor_ptr->ts_step) +
982  ".h5m");
983 
984  if (doEvalField) {
985  if constexpr (SPACE_DIM == 3) {
986  CHKERR mField.getInterface<FieldEvaluatorInterface>()
987  ->evalFEAtThePoint3D(
988  fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
989  simple->getDomainFEName(), fieldEvalData,
990  mField.get_comm_rank(), mField.get_comm_rank(), nullptr,
991  MF_EXIST, QUIET);
992  } else {
993  CHKERR mField.getInterface<FieldEvaluatorInterface>()
994  ->evalFEAtThePoint2D(
995  fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
996  simple->getDomainFEName(), fieldEvalData,
997  mField.get_comm_rank(), mField.get_comm_rank(), nullptr,
998  MF_EXIST, QUIET);
999  }
1000 
1001  if (scalarFieldPtr->size()) {
1002  auto t_temp = getFTensor0FromVec(*scalarFieldPtr);
1003  MOFEM_LOG("ThermoElasticSync", Sev::inform)
1004  << "Eval point T: " << t_temp;
1005  }
1006  if (vectorFieldPtr->size1()) {
1007  FTensor::Index<'i', SPACE_DIM> i;
1008  auto t_disp = getFTensor1FromMat<SPACE_DIM>(*vectorFieldPtr);
1009  MOFEM_LOG("ThermoElasticSync", Sev::inform)
1010  << "Eval point U magnitude: " << sqrt(t_disp(i) * t_disp(i));
1011  }
1012  if (tensorFieldPtr->size1()) {
1013  FTensor::Index<'i', SPACE_DIM> i;
1014  auto t_disp_grad =
1015  getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*tensorFieldPtr);
1016  MOFEM_LOG("ThermoElasticSync", Sev::inform)
1017  << "Eval point U_GRAD trace: " << t_disp_grad(i, i);
1018  }
1019 
1020  MOFEM_LOG_SYNCHRONISE(mField.get_comm());
1021  }
1022 
1023  MoFEMFunctionReturn(0);
1024  };
1025  auto null = boost::shared_ptr<FEMethod>();
1026  CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(), null,
1027  monitor_ptr, null);
1028  MoFEMFunctionReturn(0);
1029  };
1030 
1031  auto set_fieldsplit_preconditioner = [&](auto solver) {
1032  MoFEMFunctionBeginHot;
1033 
1034  SNES snes;
1035  CHKERR TSGetSNES(solver, &snes);
1036  KSP ksp;
1037  CHKERR SNESGetKSP(snes, &ksp);
1038  PC pc;
1039  CHKERR KSPGetPC(ksp, &pc);
1040  PetscBool is_pcfs = PETSC_FALSE;
1041  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1042 
1043  // Setup fieldsplit (block) solver - optional: yes/no
1044  if (is_pcfs == PETSC_TRUE) {
1045  auto bc_mng = mField.getInterface<BcManager>();
1046  auto is_mng = mField.getInterface<ISManager>();
1047  auto name_prb = simple->getProblemName();
1048 
1049  SmartPetscObj<IS> is_u;
1050  CHKERR is_mng->isCreateProblemFieldAndRank(name_prb, ROW, "U", 0,
1051  SPACE_DIM, is_u);
1052  SmartPetscObj<IS> is_flux;
1053  CHKERR is_mng->isCreateProblemFieldAndRank(name_prb, ROW, "FLUX", 0, 0,
1054  is_flux);
1055  SmartPetscObj<IS> is_T;
1056  CHKERR is_mng->isCreateProblemFieldAndRank(name_prb, ROW, "T", 0, 0,
1057  is_T);
1058  IS is_tmp;
1059  CHKERR ISExpand(is_T, is_flux, &is_tmp);
1060  auto is_TFlux = SmartPetscObj<IS>(is_tmp);
1061 
1062  auto is_all_bc = bc_mng->getBlockIS(name_prb, "HEATFLUX", "FLUX", 0, 0);
1063  int is_all_bc_size;
1064  CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
1065  MOFEM_LOG("ThermoElastic", Sev::inform)
1066  << "Field split block size " << is_all_bc_size;
1067  if (is_all_bc_size) {
1068  IS is_tmp2;
1069  CHKERR ISDifference(is_TFlux, is_all_bc, &is_tmp2);
1070  is_TFlux = SmartPetscObj<IS>(is_tmp2);
1071  CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
1072  is_all_bc); // boundary block
1073  }
1074 
1075  CHKERR ISSort(is_u);
1076  CHKERR ISSort(is_TFlux);
1077  CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_TFlux);
1078  CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_u);
1079  }
1080 
1081  MoFEMFunctionReturnHot(0);
1082  };
1083 
1084  auto D = createDMVector(dm);
1085  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_FORWARD);
1086  CHKERR TSSetSolution(solver, D);
1087  CHKERR TSSetFromOptions(solver);
1088 
1089  CHKERR set_section_monitor(solver);
1090  CHKERR set_fieldsplit_preconditioner(solver);
1091  CHKERR set_time_monitor(dm, solver);
1092 
1093  CHKERR TSSetUp(solver);
1094  CHKERR TSSolve(solver, NULL);
1095 
1096  MoFEMFunctionReturn(0);
1097 }
1098 //! [Solve]
1099 
1100 static char help[] = "...\n\n";
1101 
1102 int main(int argc, char *argv[]) {
1103 
1104  // Initialisation of MoFEM/PETSc and MOAB data structures
1105  const char param_file[] = "param_file.petsc";
1106  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
1107 
1108  // Add logging channel for example
1109  auto core_log = logging::core::get();
1110  core_log->add_sink(
1111  LogManager::createSink(LogManager::getStrmWorld(), "ThermoElastic"));
1112  LogManager::setLog("ThermoElastic");
1113  MOFEM_LOG_TAG("ThermoElastic", "ThermoElastic");
1114 
1115  core_log->add_sink(
1116  LogManager::createSink(LogManager::getStrmSync(), "ThermoElasticSync"));
1117  LogManager::setLog("ThermoElasticSync");
1118  MOFEM_LOG_TAG("ThermoElasticSync", "ThermoElasticSync");
1119 
1120  try {
1121 
1122  //! [Register MoFEM discrete manager in PETSc]
1123  DMType dm_name = "DMMOFEM";
1124  CHKERR DMRegister_MoFEM(dm_name);
1125  //! [Register MoFEM discrete manager in PETSc
1126 
1127  //! [Create MoAB]
1128  moab::Core mb_instance; ///< mesh database
1129  moab::Interface &moab = mb_instance; ///< mesh database interface
1130  //! [Create MoAB]
1131 
1132  //! [Create MoFEM]
1133  MoFEM::Core core(moab); ///< finite element database
1134  MoFEM::Interface &m_field = core; ///< finite element database interface
1135  //! [Create MoFEM]
1136 
1137  //! [Load mesh]
1138  Simple *simple = m_field.getInterface<Simple>();
1139  CHKERR simple->getOptions();
1140  CHKERR simple->loadFile();
1141  //! [Load mesh]
1142 
1143  //! [ThermoElasticProblem]
1144  ThermoElasticProblem ex(m_field);
1145  CHKERR ex.runProblem();
1146  //! [ThermoElasticProblem]
1147  }
1148  CATCH_ERRORS;
1149 
1150  CHKERR MoFEM::Core::Finalize();
1151 }
NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEM::NaturalBC::Assembly::LinearForm
Definition: Natural.hpp:67
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
default_coeff_expansion
double default_coeff_expansion
Definition: thermo_elastic.cpp:127
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
SIDESET
@ SIDESET
Definition: definitions.h:160
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
ThermoElasticProblem::tsSolve
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
Definition: thermo_elastic.cpp:887
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
MOFEM_LOG_SEVERITY_SYNC
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
Definition: LogManager.hpp:352
OpHdivHdiv
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, SPACE_DIM > OpHdivHdiv
[Linear elastic problem]
Definition: thermo_elastic.cpp:47
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
young_modulus
double young_modulus
Young modulus.
Definition: plastic.cpp:121
ThermoElasticProblem::ThermoElasticProblem
ThermoElasticProblem(MoFEM::Interface &m_field)
Definition: thermo_elastic.cpp:155
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
ThermoElasticProblem::scalarFieldPtr
boost::shared_ptr< VectorDouble > scalarFieldPtr
Definition: thermo_elastic.cpp:166
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
MoFEM::OpFluxRhsImpl
Definition: Natural.hpp:39
ThermoElasticProblem::BlockedParameters::getHeatCapacityPtr
auto getHeatCapacityPtr()
Definition: thermo_elastic.cpp:195
MoFEM::PipelineManager::ElementsAndOpsByDim
Definition: PipelineManager.hpp:38
MoFEM::NaturalBC::Assembly::BiLinearForm
Definition: Natural.hpp:74
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
OpHDivTemp
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 1, SPACE_DIM > OpHDivTemp
Integrate Rhs div flux base times temperature (T)
Definition: thermo_elastic.cpp:76
MoFEM::EssentialBC
Essential boundary conditions.
Definition: Essential.hpp:111
ThermoElasticOps::OpKCauchyThermoElasticity
Definition: ThermoElasticOps.hpp:9
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
ThermoElasticProblem::BlockedParameters::getHeatConductivityPtr
auto getHeatConductivityPtr()
Definition: thermo_elastic.cpp:191
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
ref_temp
double ref_temp
Definition: thermo_elastic.cpp:124
ThermoElasticProblem::fieldEvalCoords
std::array< double, SPACE_DIM > fieldEvalCoords
Definition: thermo_elastic.cpp:163
ThermoElasticProblem::bC
MoFEMErrorCode bC()
[Create common data]
Definition: thermo_elastic.cpp:556
MoFEM::OpCalculateScalarFieldValuesDot
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
Definition: UserDataOperators.hpp:273
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
ThermoElasticProblem::setupProblem
MoFEMErrorCode setupProblem()
add fields
Definition: thermo_elastic.cpp:448
MoFEM.hpp
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
bulk_modulus_K
double bulk_modulus_K
Definition: dynamic_first_order_con_law.cpp:96
ThermoElasticProblem::tensorFieldPtr
boost::shared_ptr< MatrixDouble > tensorFieldPtr
Definition: thermo_elastic.cpp:168
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
MoFEM::DisplacementCubitBcData
Definition of the displacement bc data structure.
Definition: BCData.hpp:76
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
default_poisson_ratio
double default_poisson_ratio
Definition: thermo_elastic.cpp:123
OpBase
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
MoFEM::LogManager::createSink
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
OpBaseDotT
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBaseDotT
Integrate Rhs base of temperature time heat capacity times heat rate (T)
Definition: thermo_elastic.cpp:84
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
sdf.r
int r
Definition: sdf.py:8
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::OpFluxLhsImpl
Definition: Natural.hpp:43
ROW
@ ROW
Definition: definitions.h:136
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:667
NODESET
@ NODESET
Definition: definitions.h:159
ThermoElasticOps::OpStressThermal
Definition: ThermoElasticOps.hpp:23
ThermoElasticOps.hpp
SPACE_DIM
constexpr int SPACE_DIM
Definition: thermo_elastic.cpp:17
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:178
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
MoFEM::EssentialPreProc< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:91
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
default_heat_conductivity
double default_heat_conductivity
Definition: thermo_elastic.cpp:128
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition: seepage.cpp:57
BoundaryEleOp
ThermoElasticProblem
Definition: thermo_elastic.cpp:153
MoFEM::OpCalculateHVecVectorField
Get vector field for H-div approximation.
Definition: UserDataOperators.hpp:2115
MoFEM::get_temp_meshset_ptr
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition: Templates.hpp:1886
MoFEM::FieldEvaluatorInterface
Field evaluator interface.
Definition: FieldEvaluator.hpp:21
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
ThermoElasticProblem::BlockedParameters::getDPtr
auto getDPtr()
Definition: thermo_elastic.cpp:183
double
MoFEM::OpEssentialLhsImpl
Enforce essential constrains on lhs.
Definition: Essential.hpp:81
convert.type
type
Definition: convert.py:64
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:317
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MOFEM_LOG_SYNCHRONISE
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
MoFEM::LogManager::getStrmSync
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.
Definition: LogManager.cpp:348
OpBaseDivFlux
OpBaseDotT OpBaseDivFlux
Integrate Rhs base of temperature times divergence of flux (T)
Definition: thermo_elastic.cpp:90
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition: plastic.cpp:122
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
EntData
EntitiesFieldData::EntData EntData
Definition: thermo_elastic.cpp:20
OpKCauchy
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Linear elastic problem]
Definition: thermo_elastic.cpp:34
ThermoElasticProblem::BlockedParameters::heatCapacity
double heatCapacity
Definition: thermo_elastic.cpp:181
MoFEM::AddFluxToLhsPipelineImpl
Definition: Natural.hpp:49
doEvalField
PetscBool doEvalField
Definition: incompressible_elasticity.cpp:41
ThermoElasticProblem::createCommonData
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: thermo_elastic.cpp:509
init_temp
double init_temp
Definition: thermo_elastic.cpp:125
MoFEM::LogManager::SeverityLevel
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
MoFEM::AddFluxToRhsPipelineImpl
Definition: Natural.hpp:46
ThermoElasticOps
Definition: ThermoElasticOps.hpp:7
size_symm
constexpr auto size_symm
Definition: plastic.cpp:42
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
BiLinearForm
ThermoElasticProblem::BlockedParameters::coeffExpansion
double coeffExpansion
Definition: thermo_elastic.cpp:179
OpInternalForceCauchy
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: thermo_elastic.cpp:38
ThermoElasticProblem::OPs
MoFEMErrorCode OPs()
[Boundary condition]
Definition: thermo_elastic.cpp:667
FTensor::Index< 'i', SPACE_DIM >
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:413
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1535
convert.n
n
Definition: convert.py:82
ThermoElasticProblem::BlockedParameters::getCoeffExpansionPtr
auto getCoeffExpansionPtr()
Definition: thermo_elastic.cpp:187
integration_rule
auto integration_rule
Definition: free_surface.cpp:185
ThermoElasticProblem::addMatBlockOps
MoFEMErrorCode addMatBlockOps(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, std::string block_elastic_name, std::string block_thermal_name, boost::shared_ptr< BlockedParameters > blockedParamsPtr, Sev sev)
Definition: thermo_elastic.cpp:206
MoFEM::OpUnSetBc
Definition: FormsIntegrators.hpp:49
Range
DomainEleOp
main
int main(int argc, char *argv[])
Definition: thermo_elastic.cpp:1102
MoFEM::PetscOptionsGetRealArray
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Definition: DeprecatedPetsc.hpp:192
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
save_range
auto save_range
Definition: thermo_elastic.cpp:144
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::NaturalBC::Assembly
Assembly methods.
Definition: Natural.hpp:65
MoFEM::OpEssentialRhsImpl
Enforce essential constrains on rhs.
Definition: Essential.hpp:65
ThermoElasticProblem::vectorFieldPtr
boost::shared_ptr< MatrixDouble > vectorFieldPtr
Definition: thermo_elastic.cpp:167
shear_modulus_G
double shear_modulus_G
Definition: dynamic_first_order_con_law.cpp:97
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::OpSymmetrizeTensor
Definition: UserDataOperators.hpp:1948
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
default_young_modulus
double default_young_modulus
[Essential boundary conditions (Least square approach)]
Definition: thermo_elastic.cpp:122
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
approx_order
int approx_order
Definition: test_broken_space.cpp:50
EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition: thermo_elastic.cpp:10
MoFEM::DMMoFEMTSSetMonitor
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: DMMoFEM.cpp:1056
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::OpSetBc
Set indices on entities on finite element.
Definition: FormsIntegrators.hpp:38
ThermoElasticProblem::BlockedParameters
Definition: thermo_elastic.cpp:176
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
TEMPERATURESET
@ TEMPERATURESET
Definition: definitions.h:168
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::BlockData
Definition: MeshsetsManager.cpp:755
MoFEM::OpCalculateHdivVectorDivergence
Calculate divergence of vector field.
Definition: UserDataOperators.hpp:2212
MoFEM::MoFEMSNESMonitorFields
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
Definition: SnesCtx.cpp:232
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
MoFEM::EssentialBC::Assembly
Assembly methods.
Definition: Essential.hpp:119
order
int order
Definition: thermo_elastic.cpp:134
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
OpCapacity
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpCapacity
Integrate Lhs base of temperature times (heat capacity) times base of temperature (T x T)
Definition: thermo_elastic.cpp:63
MoFEM::HeatFluxCubitBcData
Definition of the heat flux bc data structure.
Definition: BCData.hpp:427
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
default_heat_capacity
double default_heat_capacity
Definition: thermo_elastic.cpp:131
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
MoFEM::getDMSnesCtx
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition: DMMoFEM.hpp:1127
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
HEATFLUXSET
@ HEATFLUXSET
Definition: definitions.h:169
QUIET
@ QUIET
Definition: definitions.h:221
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
ThermoElasticProblem::BlockedParameters::mD
MatrixDouble mD
Definition: thermo_elastic.cpp:178
ThermoElasticProblem::mField
MoFEM::Interface & mField
Definition: thermo_elastic.cpp:160
DomainEleOp
DomainEle::UserDataOperator DomainEleOp
Definition: thermo_elastic.cpp:22
MoFEM::MeshsetsManager::getCubitMeshsetPtr
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
Definition: MeshsetsManager.cpp:578
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
MoFEM::SmartPetscObj< IS >
ThermoElasticProblem::BlockedParameters::heatConductivity
double heatConductivity
Definition: thermo_elastic.cpp:180
ThermoElasticProblem::fieldEvalData
boost::shared_ptr< FieldEvaluatorInterface::SetPtsData > fieldEvalData
Definition: thermo_elastic.cpp:164
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
OpHdivT
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,...
Definition: thermo_elastic.cpp:55
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
MF_EXIST
@ MF_EXIST
Definition: definitions.h:113
convert.int
int
Definition: convert.py:64
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
ThermoElasticProblem::runProblem
MoFEMErrorCode runProblem()
[Run problem]
Definition: thermo_elastic.cpp:436
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
help
static char help[]
[Solve]
Definition: thermo_elastic.cpp:1100
ThermoElasticProblem::doEvalField
PetscBool doEvalField
Definition: thermo_elastic.cpp:162
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698
OpHdivFlux
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 3, SPACE_DIM, 1 > OpHdivFlux
Integrating Rhs flux base (1/k) flux (FLUX)
Definition: thermo_elastic.cpp:69
PlasticOps::addMatBlockOps
MoFEMErrorCode addMatBlockOps(MoFEM::Interface &m_field, std::string block_name, Pip &pip, boost::shared_ptr< MatrixDouble > mat_D_Ptr, boost::shared_ptr< CommonData::BlockParams > mat_params_ptr, double scale_value, Sev sev)
Definition: PlasticOps.hpp:248
Generated by Doxygen 1.8.17 and hosted at University of Glasgow