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 using namespace MoFEM;
15 
16 #include <ThermalConvection.hpp>
17 #include <ThermalRadiation.hpp>
18 
19 constexpr int SPACE_DIM =
20  EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
21 
22 using EntData = EntitiesFieldData::EntData;
23 using DomainEle = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::DomainEle;
24 using DomainEleOp = DomainEle::UserDataOperator;
25 using BoundaryEle =
26  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
27 using BoundaryEleOp = BoundaryEle::UserDataOperator;
28 using PostProcEle = PostProcBrokenMeshInMoab<DomainEle>;
29 using SkinPostProcEle = PostProcBrokenMeshInMoab<BoundaryEle>;
30 using SideEle = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
31 
32 using AssemblyDomainEleOp =
33  FormsIntegrators<DomainEleOp>::Assembly<PETSC>::OpBase;
34 
35 //! [Linear elastic problem]
36 using OpKCauchy = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
37  GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM,
38  0>; //< Elastic stiffness matrix
39 using OpInternalForceCauchy =
40  FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
41  GAUSS>::OpGradTimesSymTensor<1, SPACE_DIM,
42  SPACE_DIM>; //< Elastic internal forces
43 //! [Linear elastic problem]
44 
45 //! [Thermal problem]
46 /**
47  * @brief Integrate Lhs base of flux (1/k) base of flux (FLUX x FLUX)
48  *
49  */
50 using OpHdivHdiv = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
51  GAUSS>::OpMass<3, SPACE_DIM>;
52 
53 /**
54  * @brief Integrate Lhs div of base of flux time base of temperature (FLUX x T)
55  * and transpose of it, i.e. (T x FLAX)
56  *
57  */
58 using OpHdivT = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
59  GAUSS>::OpMixDivTimesScalar<SPACE_DIM>;
60 
61 /**
62  * @brief Integrate Lhs base of temperature times (heat capacity) times base of
63  * temperature (T x T)
64  *
65  */
66 using OpCapacity = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
67  GAUSS>::OpMass<1, 1>;
68 
69 /**
70  * @brief Integrating Rhs flux base (1/k) flux (FLUX)
71  */
72 using OpHdivFlux = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
73  GAUSS>::OpBaseTimesVector<3, SPACE_DIM, 1>;
74 
75 /**
76  * @brief Integrate Rhs div flux base times temperature (T)
77  *
78  */
79 using OpHDivTemp = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
80  GAUSS>::OpMixDivTimesU<3, 1, SPACE_DIM>;
81 
82 /**
83  * @brief Integrate Rhs base of temperature time heat capacity times heat rate
84  * (T)
85  *
86  */
87 using OpBaseDotT = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
88  GAUSS>::OpBaseTimesScalar<1>;
89 
90 /**
91  * @brief Integrate Rhs base of temperature times divergence of flux (T)
92  *
93  */
94 using OpBaseDivFlux = OpBaseDotT;
95 
96 //! [Thermal problem]
97 
98 //! [Body and heat source]
99 using DomainNaturalBCRhs =
100  NaturalBC<DomainEleOp>::Assembly<PETSC>::LinearForm<GAUSS>;
101 using OpBodyForce =
102  DomainNaturalBCRhs::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, SPACE_DIM>;
103 using OpHeatSource =
104  DomainNaturalBCRhs::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 1>;
105 using DomainNaturalBCLhs =
106  NaturalBC<DomainEleOp>::Assembly<PETSC>::BiLinearForm<GAUSS>;
107 //! [Body and heat source]
108 
109 //! [Natural boundary conditions]
110 using BoundaryNaturalBC =
111  NaturalBC<BoundaryEleOp>::Assembly<PETSC>::LinearForm<GAUSS>;
112 using OpForce = BoundaryNaturalBC::OpFlux<NaturalForceMeshsets, 1, SPACE_DIM>;
113 
114 //! [Natural boundary conditions]
115 
116 //! [Essential boundary conditions (Least square approach)]
117 using OpEssentialFluxRhs =
118  EssentialBC<BoundaryEleOp>::Assembly<PETSC>::LinearForm<
119  GAUSS>::OpEssentialRhs<HeatFluxCubitBcData, 3, SPACE_DIM>;
120 using OpEssentialFluxLhs =
121  EssentialBC<BoundaryEleOp>::Assembly<PETSC>::BiLinearForm<
122  GAUSS>::OpEssentialLhs<HeatFluxCubitBcData, 3, SPACE_DIM>;
123 //! [Essential boundary conditions (Least square approach)]
124 
125 double default_young_modulus = 1;
126 double default_poisson_ratio = 0.25;
127 double ref_temp = 0.0;
128 double init_temp = 0.0;
129 
130 PetscBool is_plane_strain = PETSC_FALSE;
131 
132 double default_coeff_expansion = 1;
133 double default_heat_conductivity =
134  1; // Force / (time temperature ) or Power /
135  // (length temperature) // Time unit is hour. force unit kN
136 double default_heat_capacity = 1; // length^2/(time^2 temperature) // length is
137  // millimeter time is hour
138 
139 int order_temp = 2; //< default approximation order for the temperature field
140 int order_flux = 3; //< default approximation order for heat flux field
141 int order_disp = 3; //< default approximation order for the displacement field
142 
143 int atom_test = 0;
144 
145 int save_every = 1; //< Save every n-th step
146 PetscBool do_output_domain;
147 PetscBool do_output_skin;
148 
149 #include <ThermoElasticOps.hpp> //< additional coupling operators
150 using namespace ThermoElasticOps; //< name space of coupling operators
151 
152 using OpSetTemperatureRhs =
153  DomainNaturalBCRhs::OpFlux<SetTargetTemperature, 1, 1>;
154 using OpSetTemperatureLhs =
155  DomainNaturalBCLhs::OpFlux<SetTargetTemperature, 1, 1>;
156 
157 auto save_range = [](moab::Interface &moab, const std::string name,
158  const Range r) {
159  MoFEMFunctionBegin;
160  auto out_meshset = get_temp_meshset_ptr(moab);
161  CHKERR moab.add_entities(*out_meshset, r);
162  CHKERR moab.write_file(name.c_str(), "VTK", "", out_meshset->get_ptr(), 1);
163  MoFEMFunctionReturn(0);
164 };
165 
166 struct ThermoElasticProblem {
167 
168  ThermoElasticProblem(MoFEM::Interface &m_field) : mField(m_field) {}
169 
170  MoFEMErrorCode runProblem();
171 
172 private:
173  MoFEM::Interface &mField;
174 
175  PetscBool doEvalField;
176  std::array<double, SPACE_DIM> fieldEvalCoords;
177  boost::shared_ptr<FieldEvaluatorInterface::SetPtsData> fieldEvalData;
178 
179  boost::shared_ptr<VectorDouble> tempFieldPtr;
180  boost::shared_ptr<MatrixDouble> fluxFieldPtr;
181  boost::shared_ptr<MatrixDouble> dispFieldPtr;
182  boost::shared_ptr<MatrixDouble> dispGradPtr;
183  boost::shared_ptr<MatrixDouble> strainFieldPtr;
184  boost::shared_ptr<MatrixDouble> stressFieldPtr;
185 
186  MoFEMErrorCode setupProblem(); ///< add fields
187  MoFEMErrorCode
188  getCommandLineParameters(); //< read parameters from command line
189  MoFEMErrorCode bC(); //< read boundary conditions
190  MoFEMErrorCode OPs(); //< add operators to pipeline
191  MoFEMErrorCode tsSolve(); //< time solver
192 
193  struct BlockedParameters
194  : public boost::enable_shared_from_this<BlockedParameters> {
195  MatrixDouble mD;
196  VectorDouble coeffExpansion;
197  double heatConductivity;
198  double heatCapacity;
199 
200  inline auto getDPtr() {
201  return boost::shared_ptr<MatrixDouble>(shared_from_this(), &mD);
202  }
203 
204  inline auto getCoeffExpansionPtr() {
205  return boost::shared_ptr<VectorDouble>(shared_from_this(),
206  &coeffExpansion);
207  }
208 
209  inline auto getHeatConductivityPtr() {
210  return boost::shared_ptr<double>(shared_from_this(), &heatConductivity);
211  }
212 
213  inline auto getHeatCapacityPtr() {
214  return boost::shared_ptr<double>(shared_from_this(), &heatCapacity);
215  }
216  };
217 
218  MoFEMErrorCode addMatBlockOps(
219  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
220  std::string block_elastic_name, std::string block_thermal_name,
221  boost::shared_ptr<BlockedParameters> blockedParamsPtr, Sev sev);
222 };
223 
224 MoFEMErrorCode ThermoElasticProblem::addMatBlockOps(
225  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
226  std::string block_elastic_name, std::string block_thermal_name,
227  boost::shared_ptr<BlockedParameters> blockedParamsPtr, Sev sev) {
228  MoFEMFunctionBegin;
229 
230  struct OpMatElasticBlocks : public DomainEleOp {
231  OpMatElasticBlocks(boost::shared_ptr<MatrixDouble> m, double bulk_modulus_K,
232  double shear_modulus_G, MoFEM::Interface &m_field,
233  Sev sev,
234  std::vector<const CubitMeshSets *> meshset_vec_ptr)
235  : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), matDPtr(m),
236  bulkModulusKDefault(bulk_modulus_K),
237  shearModulusGDefault(shear_modulus_G) {
238  CHK_THROW_MESSAGE(extractElasticBlockData(m_field, meshset_vec_ptr, sev),
239  "Can not get data from block");
240  }
241 
242  MoFEMErrorCode doWork(int side, EntityType type,
243  EntitiesFieldData::EntData &data) {
244  MoFEMFunctionBegin;
245 
246  for (auto &b : blockData) {
247 
248  if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
249  CHKERR getMatDPtr(matDPtr, b.bulkModulusK, b.shearModulusG);
250  MoFEMFunctionReturnHot(0);
251  }
252  }
253 
254  CHKERR getMatDPtr(matDPtr, bulkModulusKDefault, shearModulusGDefault);
255  MoFEMFunctionReturn(0);
256  }
257 
258  private:
259  boost::shared_ptr<MatrixDouble> matDPtr;
260 
261  struct BlockData {
262  double bulkModulusK;
263  double shearModulusG;
264  Range blockEnts;
265  };
266 
267  double bulkModulusKDefault;
268  double shearModulusGDefault;
269  std::vector<BlockData> blockData;
270 
271  MoFEMErrorCode
272  extractElasticBlockData(MoFEM::Interface &m_field,
273  std::vector<const CubitMeshSets *> meshset_vec_ptr,
274  Sev sev) {
275  MoFEMFunctionBegin;
276 
277  for (auto m : meshset_vec_ptr) {
278  MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Elastic Block") << *m;
279  std::vector<double> block_data;
280  CHKERR m->getAttributes(block_data);
281  if (block_data.size() < 2) {
282  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
283  "Expected that block has two attributes");
284  }
285  auto get_block_ents = [&]() {
286  Range ents;
287  CHKERR
288  m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
289  return ents;
290  };
291 
292  double young_modulus = block_data[0];
293  double poisson_ratio = block_data[1];
294  double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
295  double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
296 
297  MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Elastic Block")
298  << m->getName() << ": E = " << young_modulus
299  << " nu = " << poisson_ratio;
300 
301  blockData.push_back(
302  {bulk_modulus_K, shear_modulus_G, get_block_ents()});
303  }
304  MOFEM_LOG_CHANNEL("WORLD");
305  MoFEMFunctionReturn(0);
306  }
307 
308  MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
309  double bulk_modulus_K, double shear_modulus_G) {
310  MoFEMFunctionBegin;
311  //! [Calculate elasticity tensor]
312  auto set_material_stiffness = [&]() {
313  FTensor::Index<'i', SPACE_DIM> i;
314  FTensor::Index<'j', SPACE_DIM> j;
315  FTensor::Index<'k', SPACE_DIM> k;
316  FTensor::Index<'l', SPACE_DIM> l;
317  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
318  double A = 1;
319  if (SPACE_DIM == 2 && !is_plane_strain) {
320  A = 2 * shear_modulus_G /
321  (bulk_modulus_K + (4. / 3.) * shear_modulus_G);
322  }
323  auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
324  t_D(i, j, k, l) =
325  2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
326  A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) * t_kd(i, j) *
327  t_kd(k, l);
328  };
329  //! [Calculate elasticity tensor]
330  constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
331  mat_D_ptr->resize(size_symm * size_symm, 1);
332  set_material_stiffness();
333  MoFEMFunctionReturn(0);
334  }
335  };
336 
337  double default_bulk_modulus_K =
338  default_young_modulus / (3 * (1 - 2 * default_poisson_ratio));
339  double default_shear_modulus_G =
340  default_young_modulus / (2 * (1 + default_poisson_ratio));
341 
342  pipeline.push_back(new OpMatElasticBlocks(
343  blockedParamsPtr->getDPtr(), default_bulk_modulus_K,
344  default_shear_modulus_G, mField, sev,
345 
346  // Get blockset using regular expression
347  mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
348 
349  (boost::format("%s(.*)") % block_elastic_name).str()
350 
351  ))
352 
353  ));
354 
355  struct OpMatThermalBlocks : public DomainEleOp {
356  OpMatThermalBlocks(boost::shared_ptr<VectorDouble> expansion_ptr,
357  boost::shared_ptr<double> conductivity_ptr,
358  boost::shared_ptr<double> capacity_ptr,
359  MoFEM::Interface &m_field, Sev sev,
360  std::vector<const CubitMeshSets *> meshset_vec_ptr)
361  : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE),
362  expansionPtr(expansion_ptr), conductivityPtr(conductivity_ptr),
363  capacityPtr(capacity_ptr) {
364  CHK_THROW_MESSAGE(extractThermallockData(m_field, meshset_vec_ptr, sev),
365  "Can not get data from block");
366  }
367 
368  MoFEMErrorCode doWork(int side, EntityType type,
369  EntitiesFieldData::EntData &data) {
370  MoFEMFunctionBegin;
371 
372  for (auto &b : blockData) {
373 
374  if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
375  *expansionPtr = b.expansion;
376  *conductivityPtr = b.conductivity;
377  *capacityPtr = b.capacity;
378  MoFEMFunctionReturnHot(0);
379  }
380  }
381 
382  *expansionPtr = VectorDouble(SPACE_DIM, default_coeff_expansion);
383  *conductivityPtr = default_heat_conductivity;
384  *capacityPtr = default_heat_capacity;
385 
386  MoFEMFunctionReturn(0);
387  }
388 
389  private:
390  struct BlockData {
391  double conductivity;
392  double capacity;
393  VectorDouble expansion;
394  Range blockEnts;
395  };
396 
397  std::vector<BlockData> blockData;
398 
399  MoFEMErrorCode
400  extractThermallockData(MoFEM::Interface &m_field,
401  std::vector<const CubitMeshSets *> meshset_vec_ptr,
402  Sev sev) {
403  MoFEMFunctionBegin;
404 
405  for (auto m : meshset_vec_ptr) {
406  MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Thermal Block") << *m;
407  std::vector<double> block_data;
408  CHKERR m->getAttributes(block_data);
409  if (block_data.size() < 3) {
410  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
411  "Expected that block has at least three attributes");
412  }
413  auto get_block_ents = [&]() {
414  Range ents;
415  CHKERR
416  m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
417  return ents;
418  };
419 
420  auto get_expansion = [&]() {
421  VectorDouble expansion(SPACE_DIM, block_data[2]);
422  if (block_data.size() > 3) {
423  expansion[1] = block_data[3];
424  }
425  if (SPACE_DIM == 3 && block_data.size() > 4) {
426  expansion[2] = block_data[4];
427  }
428  return expansion;
429  };
430 
431  auto coeff_exp_vec = get_expansion();
432 
433  MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Thermal Block")
434  << m->getName() << ": conductivity = " << block_data[0]
435  << " capacity = " << block_data[1]
436  << " expansion = " << coeff_exp_vec;
437 
438  blockData.push_back(
439  {block_data[0], block_data[1], coeff_exp_vec, get_block_ents()});
440  }
441  MOFEM_LOG_CHANNEL("WORLD");
442  MoFEMFunctionReturn(0);
443  }
444 
445  boost::shared_ptr<VectorDouble> expansionPtr;
446  boost::shared_ptr<double> conductivityPtr;
447  boost::shared_ptr<double> capacityPtr;
448  };
449 
450  pipeline.push_back(new OpMatThermalBlocks(
451  blockedParamsPtr->getCoeffExpansionPtr(),
452  blockedParamsPtr->getHeatConductivityPtr(),
453  blockedParamsPtr->getHeatCapacityPtr(), mField, sev,
454 
455  // Get blockset using regular expression
456  mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
457 
458  (boost::format("%s(.*)") % block_thermal_name).str()
459 
460  ))
461 
462  ));
463 
464  MoFEMFunctionReturn(0);
465 }
466 
467 //! [Run problem]
468 MoFEMErrorCode ThermoElasticProblem::runProblem() {
469  MoFEMFunctionBegin;
470  CHKERR getCommandLineParameters();
471  CHKERR setupProblem();
472  CHKERR bC();
473  CHKERR OPs();
474  CHKERR tsSolve();
475  MoFEMFunctionReturn(0);
476 }
477 //! [Run problem]
478 
479 //! [Get command line parameters]
480 MoFEMErrorCode ThermoElasticProblem::getCommandLineParameters() {
481  MoFEMFunctionBegin;
482 
483  auto get_command_line_parameters = [&]() {
484  MoFEMFunctionBegin;
485 
486  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order_temp,
487  PETSC_NULL);
488  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order_temp", &order_temp,
489  PETSC_NULL);
490  order_flux = order_temp + 1;
491  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order_flux", &order_flux,
492  PETSC_NULL);
493  order_disp = order_temp + 1;
494  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order_disp", &order_disp,
495  PETSC_NULL);
496  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-atom_test", &atom_test,
497  PETSC_NULL);
498  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-save_every", &save_every,
499  PETSC_NULL);
500 
501  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
502  &default_young_modulus, PETSC_NULL);
503  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
504  &default_poisson_ratio, PETSC_NULL);
505  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-plane_strain",
506  &is_plane_strain, PETSC_NULL);
507  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-coeff_expansion",
508  &default_coeff_expansion, PETSC_NULL);
509  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-ref_temp", &ref_temp,
510  PETSC_NULL);
511  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-init_temp", &init_temp,
512  PETSC_NULL);
513 
514  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-capacity",
515  &default_heat_capacity, PETSC_NULL);
516  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-conductivity",
517  &default_heat_conductivity, PETSC_NULL);
518 
519  if constexpr (SPACE_DIM == 2) {
520  do_output_domain = PETSC_TRUE;
521  do_output_skin = PETSC_FALSE;
522  } else {
523  do_output_domain = PETSC_FALSE;
524  do_output_skin = PETSC_TRUE;
525  }
526 
527  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-output_domain",
528  &do_output_domain, PETSC_NULL);
529  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-output_skin", &do_output_skin,
530  PETSC_NULL);
531 
532  MOFEM_LOG("ThermoElastic", Sev::inform)
533  << "Default Young's modulus " << default_young_modulus;
534  MOFEM_LOG("ThermoElastic", Sev::inform)
535  << "DefaultPoisson ratio " << default_poisson_ratio;
536  MOFEM_LOG("ThermoElastic", Sev::inform)
537  << "Default coeff of expansion " << default_coeff_expansion;
538  MOFEM_LOG("ThermoElastic", Sev::inform)
539  << "Default heat capacity " << default_heat_capacity;
540  MOFEM_LOG("ThermoElastic", Sev::inform)
541  << "Default heat conductivity " << default_heat_conductivity;
542 
543  MOFEM_LOG("ThermoElastic", Sev::inform)
544  << "Reference temperature " << ref_temp;
545  MOFEM_LOG("ThermoElastic", Sev::inform)
546  << "Initial temperature " << init_temp;
547 
548  MoFEMFunctionReturn(0);
549  };
550 
551  CHKERR get_command_line_parameters();
552 
553  MoFEMFunctionReturn(0);
554 }
555 //! [Get command line parameters]
556 
557 //! [Set up problem]
558 MoFEMErrorCode ThermoElasticProblem::setupProblem() {
559  MoFEMFunctionBegin;
560  Simple *simple = mField.getInterface<Simple>();
561  // Add field
562  constexpr FieldApproximationBase base = DEMKOWICZ_JACOBI_BASE;
563  // Mechanical fields
564  CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, SPACE_DIM);
565  CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE, SPACE_DIM);
566  // Temperature
567  constexpr auto flux_space = (SPACE_DIM == 2) ? HCURL : HDIV;
568  CHKERR simple->addDomainField("T", L2, base, 1);
569  CHKERR simple->addDomainField("FLUX", flux_space, base, 1);
570  CHKERR simple->addBoundaryField("FLUX", flux_space, base, 1);
571  CHKERR simple->addBoundaryField("TBC", L2, base, 1);
572 
573  CHKERR simple->setFieldOrder("U", order_disp);
574  CHKERR simple->setFieldOrder("FLUX", order_flux);
575  CHKERR simple->setFieldOrder("T", order_temp);
576  CHKERR simple->setFieldOrder("TBC", order_temp);
577 
578  CHKERR simple->setUp();
579 
580  int coords_dim = SPACE_DIM;
581  CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
582  fieldEvalCoords.data(), &coords_dim,
583  &doEvalField);
584 
585  tempFieldPtr = boost::make_shared<VectorDouble>();
586  fluxFieldPtr = boost::make_shared<MatrixDouble>();
587  dispFieldPtr = boost::make_shared<MatrixDouble>();
588  dispGradPtr = boost::make_shared<MatrixDouble>();
589  strainFieldPtr = boost::make_shared<MatrixDouble>();
590  stressFieldPtr = boost::make_shared<MatrixDouble>();
591 
592  if (doEvalField) {
593  fieldEvalData =
594  mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
595 
596  if constexpr (SPACE_DIM == 3) {
597  CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree3D(
598  fieldEvalData, simple->getDomainFEName());
599  } else {
600  CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree2D(
601  fieldEvalData, simple->getDomainFEName());
602  }
603 
604  fieldEvalData->setEvalPoints(fieldEvalCoords.data(), 1);
605  auto no_rule = [](int, int, int) { return -1; };
606 
607  auto field_eval_fe_ptr = fieldEvalData->feMethodPtr.lock();
608  field_eval_fe_ptr->getRuleHook = no_rule;
609 
610  auto block_params = boost::make_shared<BlockedParameters>();
611  auto mDPtr = block_params->getDPtr();
612  auto coeff_expansion_ptr = block_params->getCoeffExpansionPtr();
613 
614  CHKERR addMatBlockOps(field_eval_fe_ptr->getOpPtrVector(), "MAT_ELASTIC",
615  "MAT_THERMAL", block_params, Sev::verbose);
616 
617  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
618  field_eval_fe_ptr->getOpPtrVector(), {H1, HDIV});
619 
620  field_eval_fe_ptr->getOpPtrVector().push_back(
621  new OpCalculateScalarFieldValues("T", tempFieldPtr));
622  field_eval_fe_ptr->getOpPtrVector().push_back(
623  new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", fluxFieldPtr));
624  field_eval_fe_ptr->getOpPtrVector().push_back(
625  new OpCalculateVectorFieldValues<SPACE_DIM>("U", dispFieldPtr));
626  field_eval_fe_ptr->getOpPtrVector().push_back(
627  new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>("U",
628  dispGradPtr));
629  field_eval_fe_ptr->getOpPtrVector().push_back(
630  new OpSymmetrizeTensor<SPACE_DIM>(dispGradPtr, strainFieldPtr));
631  field_eval_fe_ptr->getOpPtrVector().push_back(
632  new OpStressThermal(strainFieldPtr, tempFieldPtr, mDPtr,
633  coeff_expansion_ptr, stressFieldPtr));
634  }
635 
636  MoFEMFunctionReturn(0);
637 }
638 //! [Set up problem]
639 
640 //! [Boundary condition]
641 MoFEMErrorCode ThermoElasticProblem::bC() {
642  MoFEMFunctionBegin;
643 
644  MOFEM_LOG("SYNC", Sev::noisy) << "bC";
645  MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::noisy);
646 
647  auto simple = mField.getInterface<Simple>();
648  auto bc_mng = mField.getInterface<BcManager>();
649 
650  auto get_skin = [&]() {
651  Range body_ents;
652  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
653  Skinner skin(&mField.get_moab());
654  Range skin_ents;
655  CHKERR skin.find_skin(0, body_ents, false, skin_ents);
656  return skin_ents;
657  };
658 
659  auto filter_flux_blocks = [&](auto skin, bool temp_bc = false) {
660  auto remove_cubit_blocks = [&](auto c) {
661  MoFEMFunctionBegin;
662  for (auto m :
663 
664  mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(c)
665 
666  ) {
667  Range ents;
668  CHKERR mField.get_moab().get_entities_by_dimension(
669  m->getMeshset(), SPACE_DIM - 1, ents, true);
670  skin = subtract(skin, ents);
671  }
672  MoFEMFunctionReturn(0);
673  };
674 
675  auto remove_named_blocks = [&](auto n) {
676  MoFEMFunctionBegin;
677  for (auto m : mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
678  std::regex(
679 
680  (boost::format("%s(.*)") % n).str()
681 
682  ))
683 
684  ) {
685  Range ents;
686  CHKERR mField.get_moab().get_entities_by_dimension(
687  m->getMeshset(), SPACE_DIM - 1, ents, true);
688  skin = subtract(skin, ents);
689  }
690  MoFEMFunctionReturn(0);
691  };
692  if (!temp_bc) {
693  CHK_THROW_MESSAGE(remove_cubit_blocks(NODESET | TEMPERATURESET),
694  "remove_cubit_blocks");
695  CHK_THROW_MESSAGE(remove_named_blocks("TEMPERATURE"),
696  "remove_named_blocks");
697  }
698  CHK_THROW_MESSAGE(remove_cubit_blocks(SIDESET | HEATFLUXSET),
699  "remove_cubit_blocks");
700  CHK_THROW_MESSAGE(remove_named_blocks("HEATFLUX"), "remove_named_blocks");
701  CHK_THROW_MESSAGE(remove_named_blocks("CONVECTION"), "remove_named_blocks");
702  CHK_THROW_MESSAGE(remove_named_blocks("RADIATION"), "remove_named_blocks");
703  return skin;
704  };
705 
706  auto filter_true_skin = [&](auto skin) {
707  Range boundary_ents;
708  ParallelComm *pcomm =
709  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
710  CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
711  PSTATUS_NOT, -1, &boundary_ents);
712  return boundary_ents;
713  };
714 
715  auto remove_flux_ents = filter_true_skin(filter_flux_blocks(get_skin()));
716  auto remove_temp_bc_ents =
717  filter_true_skin(filter_flux_blocks(get_skin(), true));
718 
719  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
720  remove_flux_ents);
721  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
722  remove_temp_bc_ents);
723 
724  MOFEM_LOG("SYNC", Sev::noisy) << remove_flux_ents << endl;
725  MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::noisy);
726 
727  MOFEM_LOG("SYNC", Sev::noisy) << remove_temp_bc_ents << endl;
728  MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::noisy);
729 
730 #ifndef NDEBUG
731 
732  CHKERR save_range(
733  mField.get_moab(),
734  (boost::format("flux_remove_%d.vtk") % mField.get_comm_rank()).str(),
735  remove_flux_ents);
736 
737  CHKERR save_range(
738  mField.get_moab(),
739  (boost::format("temp_bc_remove_%d.vtk") % mField.get_comm_rank()).str(),
740  remove_temp_bc_ents);
741 
742 #endif
743 
744  CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
745  simple->getProblemName(), "FLUX", remove_flux_ents);
746  CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
747  simple->getProblemName(), "TBC", remove_temp_bc_ents);
748 
749  auto set_init_temp = [](boost::shared_ptr<FieldEntity> field_entity_ptr) {
750  field_entity_ptr->getEntFieldData()[0] = init_temp;
751  return 0;
752  };
753  CHKERR mField.getInterface<FieldBlas>()->fieldLambdaOnEntities(set_init_temp,
754  "T");
755 
756  CHKERR bc_mng->removeBlockDOFsOnEntities<DisplacementCubitBcData>(
757  simple->getProblemName(), "U");
758  CHKERR bc_mng->pushMarkDOFsOnEntities<HeatFluxCubitBcData>(
759  simple->getProblemName(), "FLUX", false);
760 
761  MoFEMFunctionReturn(0);
762 }
763 //! [Boundary condition]
764 
765 //! [Push operators to pipeline]
766 MoFEMErrorCode ThermoElasticProblem::OPs() {
767  MoFEMFunctionBegin;
768 
769  MOFEM_LOG("SYNC", Sev::noisy) << "OPs";
770  MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::noisy);
771 
772  auto pipeline_mng = mField.getInterface<PipelineManager>();
773  auto simple = mField.getInterface<Simple>();
774  auto bc_mng = mField.getInterface<BcManager>();
775 
776  auto boundary_marker =
777  bc_mng->getMergedBlocksMarker(vector<string>{"HEATFLUX"});
778 
779  auto integration_rule = [](int, int, int approx_order) {
780  return 2 * approx_order;
781  };
782  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
783  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
784  CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
785  CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
786 
787  auto block_params = boost::make_shared<BlockedParameters>();
788  auto mDPtr = block_params->getDPtr();
789  auto coeff_expansion_ptr = block_params->getCoeffExpansionPtr();
790  auto heat_conductivity_ptr = block_params->getHeatConductivityPtr();
791  auto heat_capacity_ptr = block_params->getHeatCapacityPtr();
792 
793  // Default time scaling of BCs and sources, from command line arguments
794  auto time_scale =
795  boost::make_shared<TimeScale>("", false, [](const double) { return 1; });
796  auto def_time_scale = [time_scale](const double time) {
797  return (!time_scale->argFileScale) ? time : 1;
798  };
799  auto def_file_name = [time_scale](const std::string &&name) {
800  return (!time_scale->argFileScale) ? name : "";
801  };
802 
803  // Files which define scaling for separate variables, if provided
804  auto time_bodyforce_scaling = boost::make_shared<TimeScale>(
805  def_file_name("bodyforce_scale.txt"), false, def_time_scale);
806  auto time_heatsource_scaling = boost::make_shared<TimeScale>(
807  def_file_name("heatsource_scale.txt"), false, def_time_scale);
808  auto time_temperature_scaling = boost::make_shared<TimeScale>(
809  def_file_name("temperature_bc_scale.txt"), false, def_time_scale);
810  auto time_displacement_scaling = boost::make_shared<TimeScale>(
811  def_file_name("displacement_bc_scale.txt"), false, def_time_scale);
812  auto time_flux_scaling = boost::make_shared<TimeScale>(
813  def_file_name("flux_bc_scale.txt"), false, def_time_scale);
814  auto time_force_scaling = boost::make_shared<TimeScale>(
815  def_file_name("force_bc_scale.txt"), false, def_time_scale);
816 
817  auto add_domain_rhs_ops = [&](auto &pipeline) {
818  MoFEMFunctionBegin;
819  CHKERR addMatBlockOps(pipeline, "MAT_ELASTIC", "MAT_THERMAL", block_params,
820  Sev::inform);
821  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pipeline, {H1, HDIV});
822 
823  auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
824  auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
825  auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
826 
827  auto vec_temp_ptr = boost::make_shared<VectorDouble>();
828  auto vec_temp_dot_ptr = boost::make_shared<VectorDouble>();
829  auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
830  auto vec_temp_div_ptr = boost::make_shared<VectorDouble>();
831 
832  pipeline.push_back(new OpCalculateScalarFieldValues("T", vec_temp_ptr));
833  pipeline.push_back(
834  new OpCalculateScalarFieldValuesDot("T", vec_temp_dot_ptr));
835  pipeline.push_back(new OpCalculateHdivVectorDivergence<3, SPACE_DIM>(
836  "FLUX", vec_temp_div_ptr));
837  pipeline.push_back(
838  new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
839 
840  pipeline.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
841  "U", mat_grad_ptr));
842  pipeline.push_back(
843  new OpSymmetrizeTensor<SPACE_DIM>(mat_grad_ptr, mat_strain_ptr));
844  pipeline.push_back(new OpStressThermal(mat_strain_ptr, vec_temp_ptr, mDPtr,
845  coeff_expansion_ptr,
846  mat_stress_ptr));
847 
848  pipeline.push_back(new OpInternalForceCauchy(
849  "U", mat_stress_ptr,
850  [](double, double, double) constexpr { return 1; }));
851 
852  auto resistance = [heat_conductivity_ptr](const double, const double,
853  const double) {
854  return (1. / (*heat_conductivity_ptr));
855  };
856  // negative value is consequence of symmetric system
857  auto capacity = [heat_capacity_ptr](const double, const double,
858  const double) {
859  return -(*heat_capacity_ptr);
860  };
861  auto unity = [](const double, const double, const double) constexpr {
862  return -1.;
863  };
864  pipeline.push_back(new OpHdivFlux("FLUX", mat_flux_ptr, resistance));
865  pipeline.push_back(new OpHDivTemp("FLUX", vec_temp_ptr, unity));
866  pipeline.push_back(new OpBaseDivFlux("T", vec_temp_div_ptr, unity));
867  pipeline.push_back(new OpBaseDotT("T", vec_temp_dot_ptr, capacity));
868 
869  // Set source terms
870  CHKERR DomainNaturalBCRhs::AddFluxToPipeline<OpHeatSource>::add(
871  pipeline, mField, "T", {time_scale, time_heatsource_scaling},
872  "HEAT_SOURCE", Sev::inform);
873  CHKERR DomainNaturalBCRhs::AddFluxToPipeline<OpBodyForce>::add(
874  pipeline, mField, "U", {time_scale, time_bodyforce_scaling},
875  "BODY_FORCE", Sev::inform);
876  CHKERR DomainNaturalBCRhs::AddFluxToPipeline<OpSetTemperatureRhs>::add(
877  pipeline, mField, "T", vec_temp_ptr, "SETTEMP", Sev::inform);
878 
879  MoFEMFunctionReturn(0);
880  };
881 
882  auto add_domain_lhs_ops = [&](auto &pipeline) {
883  MoFEMFunctionBegin;
884  CHKERR addMatBlockOps(pipeline, "MAT_ELASTIC", "MAT_THERMAL", block_params,
885  Sev::verbose);
886  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pipeline, {H1, HDIV});
887 
888  pipeline.push_back(new OpKCauchy("U", "U", mDPtr));
889  pipeline.push_back(new ThermoElasticOps::OpKCauchyThermoElasticity(
890  "U", "T", mDPtr, coeff_expansion_ptr));
891 
892  auto resistance = [heat_conductivity_ptr](const double, const double,
893  const double) {
894  return (1. / (*heat_conductivity_ptr));
895  };
896  auto capacity = [heat_capacity_ptr](const double, const double,
897  const double) {
898  return -(*heat_capacity_ptr);
899  };
900  pipeline.push_back(new OpHdivHdiv("FLUX", "FLUX", resistance));
901  pipeline.push_back(
902  new OpHdivT("FLUX", "T", []() constexpr { return -1; }, true));
903 
904  auto op_capacity = new OpCapacity("T", "T", capacity);
905  op_capacity->feScalingFun = [](const FEMethod *fe_ptr) {
906  return fe_ptr->ts_a;
907  };
908  pipeline.push_back(op_capacity);
909 
910  auto vec_temp_ptr = boost::make_shared<VectorDouble>();
911  pipeline.push_back(new OpCalculateScalarFieldValues("T", vec_temp_ptr));
912  CHKERR DomainNaturalBCLhs::AddFluxToPipeline<OpSetTemperatureLhs>::add(
913  pipeline, mField, "T", vec_temp_ptr, "SETTEMP", Sev::verbose);
914 
915  MoFEMFunctionReturn(0);
916  };
917 
918  auto add_boundary_rhs_ops = [&](auto &pipeline) {
919  MoFEMFunctionBegin;
920 
921  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pipeline, {HDIV});
922 
923  CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpForce>::add(
924  pipeline, mField, "U", {time_scale, time_force_scaling}, "FORCE",
925  Sev::inform);
926 
927  // Temperature BC
928 
929  using OpTemperatureBC =
930  BoundaryNaturalBC::OpFlux<NaturalTemperatureMeshsets, 3, SPACE_DIM>;
931  CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpTemperatureBC>::add(
932  pipeline, mField, "FLUX", {time_scale, time_temperature_scaling},
933  "TEMPERATURE", Sev::inform);
934 
935  // Note: fluxes for temperature should be aggregated. Here separate is
936  // NaturalMeshsetType<HEATFLUXSET>, we should add version with BLOCKSET,
937  // convection, see example tutorials/vec-0/src/NaturalBoundaryBC.hpp.
938  // and vec-0/elastic.cpp
939 
940  using OpFluxBC =
941  BoundaryNaturalBC::OpFlux<NaturalMeshsetType<HEATFLUXSET>, 1, 1>;
942  CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpFluxBC>::add(
943  pipeline, mField, "TBC", {time_scale, time_flux_scaling}, "FLUX",
944  Sev::inform);
945 
946  using T = NaturalBC<BoundaryEleOp>::Assembly<PETSC>::LinearForm<GAUSS>;
947  using OpConvectionRhsBC =
948  T::OpFlux<ThermoElasticOps::ConvectionBcType<BLOCKSET>, 1, 1>;
949  using OpRadiationRhsBC =
950  T::OpFlux<ThermoElasticOps::RadiationBcType<BLOCKSET>, 1, 1>;
951  auto temp_bc_ptr = boost::make_shared<VectorDouble>();
952  pipeline.push_back(new OpCalculateScalarFieldValues("TBC", temp_bc_ptr));
953  T::AddFluxToPipeline<OpConvectionRhsBC>::add(
954  pipeline, mField, "TBC", temp_bc_ptr, "CONVECTION", Sev::inform);
955  T::AddFluxToPipeline<OpRadiationRhsBC>::add(
956  pipeline, mField, "TBC", temp_bc_ptr, "RADIATION", Sev::inform);
957 
958  auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
959  pipeline.push_back(
960  new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
961 
962  // This is temporary implementation. It be moved to LinearFormsIntegrators,
963  // however for hybridised case, what is goal of this changes, such function
964  // is implemented for fluxes in broken space. Thus ultimately this operator
965  // would be not needed.
966 
967  struct OpTBCTimesNormalFlux
968  : public FormsIntegrators<BoundaryEleOp>::Assembly<PETSC>::OpBase {
969 
970  using OP = FormsIntegrators<BoundaryEleOp>::Assembly<PETSC>::OpBase;
971 
972  OpTBCTimesNormalFlux(const std::string field_name,
973  boost::shared_ptr<MatrixDouble> vec,
974  boost::shared_ptr<Range> ents_ptr = nullptr)
975  : OP(field_name, field_name, OP::OPROW, ents_ptr), sourceVec(vec) {}
976 
977  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data) {
978  MoFEMFunctionBegin;
979  FTENSOR_INDEX(SPACE_DIM, i);
980  // get integration weights
981  auto t_w = OP::getFTensor0IntegrationWeight();
982  // get base function values on rows
983  auto t_row_base = row_data.getFTensor0N();
984  // get normal vector
985  auto t_normal = OP::getFTensor1NormalsAtGaussPts();
986  // get vector values
987  auto t_vec = getFTensor1FromMat<SPACE_DIM, 1>(*sourceVec);
988  // loop over integration points
989  for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
990  // take into account Jacobian
991  const double alpha = t_w * (t_vec(i) * t_normal(i));
992  // loop over rows base functions
993  int rr = 0;
994  for (; rr != OP::nbRows; ++rr) {
995  OP::locF[rr] += alpha * t_row_base;
996  ++t_row_base;
997  }
998  for (; rr < OP::nbRowBaseFunctions; ++rr)
999  ++t_row_base;
1000  ++t_w; // move to another integration weight
1001  ++t_vec;
1002  ++t_normal;
1003  }
1004  EntityType fe_type = OP::getNumeredEntFiniteElementPtr()->getEntType();
1005  if (fe_type == MBTRI) {
1006  OP::locF /= 2;
1007  }
1008  MoFEMFunctionReturn(0);
1009  }
1010 
1011  protected:
1012  boost::shared_ptr<MatrixDouble> sourceVec;
1013  };
1014  pipeline.push_back(new OpTBCTimesNormalFlux("TBC", mat_flux_ptr));
1015 
1016  struct OpBaseNormalTimesTBC
1017  : public FormsIntegrators<BoundaryEleOp>::Assembly<PETSC>::OpBase {
1018 
1019  using OP = FormsIntegrators<BoundaryEleOp>::Assembly<PETSC>::OpBase;
1020 
1021  OpBaseNormalTimesTBC(const std::string field_name,
1022  boost::shared_ptr<VectorDouble> vec,
1023  boost::shared_ptr<Range> ents_ptr = nullptr)
1024  : OP(field_name, field_name, OP::OPROW, ents_ptr), sourceVec(vec) {}
1025 
1026  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data) {
1027  MoFEMFunctionBegin;
1028  FTENSOR_INDEX(SPACE_DIM, i);
1029  // get integration weights
1030  auto t_w = OP::getFTensor0IntegrationWeight();
1031  // get base function values on rows
1032  auto t_row_base = row_data.getFTensor1N<3>();
1033  // get normal vector
1034  auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1035  // get vector values
1036  auto t_vec = getFTensor0FromVec(*sourceVec);
1037  // loop over integration points
1038  for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
1039  // take into account Jacobian
1040  const double alpha = t_w * t_vec;
1041  // loop over rows base functions
1042  int rr = 0;
1043  for (; rr != OP::nbRows; ++rr) {
1044  OP::locF[rr] += alpha * (t_row_base(i) * t_normal(i));
1045  ++t_row_base;
1046  }
1047  for (; rr < OP::nbRowBaseFunctions; ++rr)
1048  ++t_row_base;
1049  ++t_w; // move to another integration weight
1050  ++t_vec;
1051  ++t_normal;
1052  }
1053  EntityType fe_type = OP::getNumeredEntFiniteElementPtr()->getEntType();
1054  if (fe_type == MBTRI) {
1055  OP::locF /= 2;
1056  }
1057  MoFEMFunctionReturn(0);
1058  }
1059 
1060  protected:
1061  boost::shared_ptr<VectorDouble> sourceVec;
1062  };
1063 
1064  pipeline.push_back(new OpBaseNormalTimesTBC("FLUX", temp_bc_ptr));
1065 
1066  MoFEMFunctionReturn(0);
1067  };
1068 
1069  auto add_boundary_lhs_ops = [&](auto &pipeline) {
1070  MoFEMFunctionBegin;
1071 
1072  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pipeline, {HDIV});
1073 
1074  using T =
1075  NaturalBC<BoundaryEleOp>::template Assembly<PETSC>::BiLinearForm<GAUSS>;
1076  using OpConvectionLhsBC =
1077  T::OpFlux<ThermoElasticOps::ConvectionBcType<BLOCKSET>, 1, 1>;
1078  using OpRadiationLhsBC =
1079  T::OpFlux<ThermoElasticOps::RadiationBcType<BLOCKSET>, 1, 1>;
1080  auto temp_bc_ptr = boost::make_shared<VectorDouble>();
1081  pipeline.push_back(new OpCalculateScalarFieldValues("TBC", temp_bc_ptr));
1082  T::AddFluxToPipeline<OpConvectionLhsBC>::add(pipeline, mField, "TBC", "TBC",
1083  "CONVECTION", Sev::verbose);
1084  T::AddFluxToPipeline<OpRadiationLhsBC>::add(
1085  pipeline, mField, "TBC", "TBC", temp_bc_ptr, "RADIATION", Sev::verbose);
1086 
1087  // This is temporary implementation. It be moved to LinearFormsIntegrators,
1088  // however for hybridised case, what is goal of this changes, such function
1089  // is implemented for fluxes in broken space. Thus ultimately this operator
1090  // would be not needed.
1091 
1092  struct OpTBCTimesNormalFlux
1093  : public FormsIntegrators<BoundaryEleOp>::Assembly<PETSC>::OpBase {
1094 
1095  using OP = FormsIntegrators<BoundaryEleOp>::Assembly<PETSC>::OpBase;
1096 
1097  OpTBCTimesNormalFlux(const std::string row_field_name,
1098  const std::string col_field_name,
1099  boost::shared_ptr<Range> ents_ptr = nullptr)
1100  : OP(row_field_name, col_field_name, OP::OPROWCOL, ents_ptr) {
1101  this->sYmm = false;
1102  this->assembleTranspose = true;
1103  this->onlyTranspose = false;
1104  }
1105 
1106  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1107  EntitiesFieldData::EntData &col_data) {
1108  MoFEMFunctionBegin;
1109 
1110  FTENSOR_INDEX(SPACE_DIM, i);
1111 
1112  // get integration weights
1113  auto t_w = OP::getFTensor0IntegrationWeight();
1114  // get base function values on rows
1115  auto t_row_base = row_data.getFTensor0N();
1116  // get normal vector
1117  auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1118  // loop over integration points
1119  for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
1120  // loop over rows base functions
1121  auto a_mat_ptr = &*OP::locMat.data().begin();
1122  int rr = 0;
1123  for (; rr != OP::nbRows; rr++) {
1124  // take into account Jacobian
1125  const double alpha = t_w * t_row_base;
1126  // get column base functions values at gauss point gg
1127  auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
1128  // loop over columns
1129  for (int cc = 0; cc != OP::nbCols; cc++) {
1130  // calculate element of local matrix
1131  // using L2 norm of normal vector for convenient area calculation
1132  // for quads, tris etc.
1133  *a_mat_ptr += alpha * (t_col_base(i) * t_normal(i));
1134  ++t_col_base;
1135  ++a_mat_ptr;
1136  }
1137  ++t_row_base;
1138  }
1139  for (; rr < OP::nbRowBaseFunctions; ++rr)
1140  ++t_row_base;
1141  ++t_normal;
1142  ++t_w; // move to another integration weight
1143  }
1144  EntityType fe_type = OP::getNumeredEntFiniteElementPtr()->getEntType();
1145  if (fe_type == MBTRI) {
1146  OP::locMat /= 2;
1147  }
1148  MoFEMFunctionReturn(0);
1149  }
1150  };
1151 
1152  pipeline.push_back(new OpTBCTimesNormalFlux("TBC", "FLUX"));
1153 
1154  MoFEMFunctionReturn(0);
1155  };
1156 
1157  // Set BCs by eliminating degrees of freedom
1158  auto get_bc_hook_rhs = [&]() {
1159  EssentialPreProc<DisplacementCubitBcData> hook(
1160  mField, pipeline_mng->getDomainRhsFE(),
1161  {time_scale, time_displacement_scaling});
1162  return hook;
1163  };
1164  auto get_bc_hook_lhs = [&]() {
1165  EssentialPreProc<DisplacementCubitBcData> hook(
1166  mField, pipeline_mng->getDomainLhsFE(),
1167  {time_scale, time_displacement_scaling});
1168  return hook;
1169  };
1170 
1171  pipeline_mng->getDomainRhsFE()->preProcessHook = get_bc_hook_rhs();
1172  pipeline_mng->getDomainLhsFE()->preProcessHook = get_bc_hook_lhs();
1173 
1174  CHKERR add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
1175  CHKERR add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
1176  CHKERR add_boundary_rhs_ops(pipeline_mng->getOpBoundaryRhsPipeline());
1177  CHKERR add_boundary_lhs_ops(pipeline_mng->getOpBoundaryLhsPipeline());
1178 
1179  MoFEMFunctionReturn(0);
1180 }
1181 //! [Push operators to pipeline]
1182 
1183 //! [Solve]
1184 MoFEMErrorCode ThermoElasticProblem::tsSolve() {
1185  MoFEMFunctionBegin;
1186 
1187  Simple *simple = mField.getInterface<Simple>();
1188  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
1189  ISManager *is_manager = mField.getInterface<ISManager>();
1190 
1191  auto dm = simple->getDM();
1192  auto solver = pipeline_mng->createTSIM();
1193  auto snes_ctx_ptr = getDMSnesCtx(dm);
1194 
1195  auto set_section_monitor = [&](auto solver) {
1196  MoFEMFunctionBegin;
1197  SNES snes;
1198  CHKERR TSGetSNES(solver, &snes);
1199  CHKERR SNESMonitorSet(snes,
1200  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
1201  void *))MoFEMSNESMonitorFields,
1202  (void *)(snes_ctx_ptr.get()), nullptr);
1203  MoFEMFunctionReturn(0);
1204  };
1205 
1206  auto create_post_process_elements = [&]() {
1207  auto block_params = boost::make_shared<BlockedParameters>();
1208  auto mDPtr = block_params->getDPtr();
1209  auto coeff_expansion_ptr = block_params->getCoeffExpansionPtr();
1210  auto u_ptr = boost::make_shared<MatrixDouble>();
1211  auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
1212  auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
1213  auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
1214  auto vec_temp_ptr = boost::make_shared<VectorDouble>();
1215  auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
1216 
1217  auto push_domain_ops = [&](auto &pp_fe) {
1218  MoFEMFunctionBegin;
1219  auto &pip = pp_fe->getOpPtrVector();
1220 
1221  CHKERR addMatBlockOps(pip, "MAT_ELASTIC", "MAT_THERMAL", block_params,
1222  Sev::verbose);
1223 
1224  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV});
1225 
1226  pip.push_back(new OpCalculateScalarFieldValues("T", vec_temp_ptr));
1227  pip.push_back(
1228  new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
1229 
1230  pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
1231  pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
1232  "U", mat_grad_ptr));
1233  pip.push_back(
1234  new OpSymmetrizeTensor<SPACE_DIM>(mat_grad_ptr, mat_strain_ptr));
1235  pip.push_back(new OpStressThermal(mat_strain_ptr, vec_temp_ptr, mDPtr,
1236  coeff_expansion_ptr, mat_stress_ptr));
1237  MoFEMFunctionReturn(0);
1238  };
1239 
1240  auto push_post_proc_ops = [&](auto &pp_fe) {
1241  MoFEMFunctionBegin;
1242  auto &pip = pp_fe->getOpPtrVector();
1243  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
1244 
1245  pip.push_back(
1246 
1247  new OpPPMap(
1248 
1249  pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
1250 
1251  {{"T", vec_temp_ptr}},
1252 
1253  {{"U", u_ptr}, {"FLUX", mat_flux_ptr}},
1254 
1255  {},
1256 
1257  {{"STRAIN", mat_strain_ptr}, {"STRESS", mat_stress_ptr}}
1258 
1259  )
1260 
1261  );
1262  MoFEMFunctionReturn(0);
1263  };
1264 
1265  auto domain_post_proc = [&]() {
1266  if (do_output_domain == PETSC_FALSE)
1267  return boost::shared_ptr<PostProcEle>();
1268  auto pp_fe = boost::make_shared<PostProcEle>(mField);
1269  CHK_MOAB_THROW(push_domain_ops(pp_fe),
1270  "push domain ops to domain element");
1271  CHK_MOAB_THROW(push_post_proc_ops(pp_fe),
1272  "push post proc ops to domain element");
1273  return pp_fe;
1274  };
1275 
1276  auto skin_post_proc = [&]() {
1277  if (do_output_skin == PETSC_FALSE)
1278  return boost::shared_ptr<SkinPostProcEle>();
1279  auto pp_fe = boost::make_shared<SkinPostProcEle>(mField);
1280  auto simple = mField.getInterface<Simple>();
1281  auto op_side = new OpLoopSide<SideEle>(mField, simple->getDomainFEName(),
1282  SPACE_DIM, Sev::verbose);
1283  CHK_MOAB_THROW(push_domain_ops(op_side),
1284  "push domain ops to side element");
1285  pp_fe->getOpPtrVector().push_back(op_side);
1286  CHK_MOAB_THROW(push_post_proc_ops(pp_fe),
1287  "push post proc ops to skin element");
1288  return pp_fe;
1289  };
1290 
1291  return std::make_pair(domain_post_proc(), skin_post_proc());
1292  };
1293 
1294  auto monitor_ptr = boost::make_shared<FEMethod>();
1295 
1296  auto [domain_post_proc_fe, skin_post_proc_fe] =
1297  create_post_process_elements();
1298 
1299  auto set_time_monitor = [&](auto dm, auto solver) {
1300  MoFEMFunctionBegin;
1301  monitor_ptr->preProcessHook = [&]() {
1302  MoFEMFunctionBegin;
1303 
1304  if (save_every && (monitor_ptr->ts_step % save_every == 0)) {
1305  if (do_output_domain) {
1306  CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
1307  domain_post_proc_fe,
1308  monitor_ptr->getCacheWeakPtr());
1309  CHKERR domain_post_proc_fe->writeFile(
1310  "out_" + boost::lexical_cast<std::string>(monitor_ptr->ts_step) +
1311  ".h5m");
1312  }
1313  if (do_output_skin) {
1314  CHKERR DMoFEMLoopFiniteElements(dm, simple->getBoundaryFEName(),
1315  skin_post_proc_fe,
1316  monitor_ptr->getCacheWeakPtr());
1317  CHKERR skin_post_proc_fe->writeFile(
1318  "out_skin_" +
1319  boost::lexical_cast<std::string>(monitor_ptr->ts_step) + ".h5m");
1320  }
1321  }
1322 
1323  if (doEvalField) {
1324 
1325  if constexpr (SPACE_DIM == 3) {
1326  CHKERR mField.getInterface<FieldEvaluatorInterface>()
1327  ->evalFEAtThePoint3D(
1328  fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
1329  simple->getDomainFEName(), fieldEvalData,
1330  mField.get_comm_rank(), mField.get_comm_rank(), nullptr,
1331  MF_EXIST, QUIET);
1332  } else {
1333  CHKERR mField.getInterface<FieldEvaluatorInterface>()
1334  ->evalFEAtThePoint2D(
1335  fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
1336  simple->getDomainFEName(), fieldEvalData,
1337  mField.get_comm_rank(), mField.get_comm_rank(), nullptr,
1338  MF_EXIST, QUIET);
1339  }
1340 
1341  if (atom_test) {
1342  auto eval_num_vec =
1343  createVectorMPI(mField.get_comm(), PETSC_DECIDE, 1);
1344  CHKERR VecZeroEntries(eval_num_vec);
1345  if (tempFieldPtr->size()) {
1346  CHKERR VecSetValue(eval_num_vec, 0, 1, ADD_VALUES);
1347  }
1348  CHKERR VecAssemblyBegin(eval_num_vec);
1349  CHKERR VecAssemblyEnd(eval_num_vec);
1350 
1351  double eval_num;
1352  CHKERR VecSum(eval_num_vec, &eval_num);
1353  if (!(int)eval_num) {
1354  SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1355  "atom test %d failed: did not find elements to evaluate "
1356  "the field, check the coordinates",
1357  atom_test);
1358  }
1359  }
1360 
1361  if (tempFieldPtr->size()) {
1362  auto t_temp = getFTensor0FromVec(*tempFieldPtr);
1363  MOFEM_LOG("ThermoElasticSync", Sev::inform)
1364  << "Eval point T: " << t_temp;
1365  if (atom_test && fabs(monitor_ptr->ts_t - 10) < 1e-12) {
1366  if (atom_test <= 3 && fabs(t_temp - 554.48) > 1e-2) {
1367  SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1368  "atom test %d failed: wrong temperature value",
1369  atom_test);
1370  }
1371  if (atom_test == 4 && fabs(t_temp - 250) > 1e-2) {
1372  SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1373  "atom test %d failed: wrong temperature value",
1374  atom_test);
1375  }
1376  }
1377  }
1378  if (fluxFieldPtr->size1()) {
1379  FTensor::Index<'i', SPACE_DIM> i;
1380  auto t_flux = getFTensor1FromMat<SPACE_DIM>(*fluxFieldPtr);
1381  auto flux_mag = sqrt(t_flux(i) * t_flux(i));
1382  MOFEM_LOG("ThermoElasticSync", Sev::inform)
1383  << "Eval point FLUX magnitude: " << flux_mag;
1384  if (atom_test && fabs(monitor_ptr->ts_t - 10) < 1e-12) {
1385  if (atom_test <= 3 && fabs(flux_mag - 27008.0) > 2e1) {
1386  SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1387  "atom test %d failed: wrong flux value", atom_test);
1388  }
1389  if (atom_test == 4 && fabs(flux_mag - 150e3) > 1e-1) {
1390  SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1391  "atom test %d failed: wrong flux value", atom_test);
1392  }
1393  }
1394  }
1395  if (dispFieldPtr->size1()) {
1396  FTensor::Index<'i', SPACE_DIM> i;
1397  auto t_disp = getFTensor1FromMat<SPACE_DIM>(*dispFieldPtr);
1398  auto disp_mag = sqrt(t_disp(i) * t_disp(i));
1399  MOFEM_LOG("ThermoElasticSync", Sev::inform)
1400  << "Eval point U magnitude: " << disp_mag;
1401  if (atom_test && fabs(monitor_ptr->ts_t - 10) < 1e-12) {
1402  if (atom_test == 1 && fabs(disp_mag - 0.00345) > 1e-5) {
1403  SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1404  "atom test %d failed: wrong displacement value",
1405  atom_test);
1406  }
1407  if ((atom_test == 2 || atom_test == 3) &&
1408  fabs(disp_mag - 0.00265) > 1e-5) {
1409  SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1410  "atom test %d failed: wrong displacement value",
1411  atom_test);
1412  }
1413  }
1414  }
1415  if (strainFieldPtr->size1()) {
1416  FTensor::Index<'i', SPACE_DIM> i;
1417  auto t_strain =
1418  getFTensor2SymmetricFromMat<SPACE_DIM>(*strainFieldPtr);
1419  auto t_strain_trace = t_strain(i, i);
1420  if (atom_test && fabs(monitor_ptr->ts_t - 10) < 1e-12) {
1421  if (atom_test == 1 && fabs(t_strain_trace - 0.00679) > 1e-5) {
1422  SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1423  "atom test %d failed: wrong strain value", atom_test);
1424  }
1425  if ((atom_test == 2 || atom_test == 3) &&
1426  fabs(t_strain_trace - 0.00522) > 1e-5) {
1427  SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1428  "atom test %d failed: wrong strain value", atom_test);
1429  }
1430  }
1431  }
1432  if (stressFieldPtr->size1()) {
1433  auto t_stress =
1434  getFTensor2SymmetricFromMat<SPACE_DIM>(*stressFieldPtr);
1435  auto von_mises_stress = std::sqrt(
1436  0.5 *
1437  ((t_stress(0, 0) - t_stress(1, 1)) *
1438  (t_stress(0, 0) - t_stress(1, 1)) +
1439  (SPACE_DIM == 3 ? (t_stress(1, 1) - t_stress(2, 2)) *
1440  (t_stress(1, 1) - t_stress(2, 2))
1441  : 0) +
1442  (SPACE_DIM == 3 ? (t_stress(2, 2) - t_stress(0, 0)) *
1443  (t_stress(2, 2) - t_stress(0, 0))
1444  : 0) +
1445  6 * (t_stress(0, 1) * t_stress(0, 1) +
1446  (SPACE_DIM == 3 ? t_stress(1, 2) * t_stress(1, 2) : 0) +
1447  (SPACE_DIM == 3 ? t_stress(2, 0) * t_stress(2, 0) : 0))));
1448  MOFEM_LOG("ThermoElasticSync", Sev::inform)
1449  << "Eval point von Mises Stress: " << von_mises_stress;
1450  if (atom_test && fabs(monitor_ptr->ts_t - 10) < 1e-12) {
1451  if (atom_test == 1 && fabs(von_mises_stress - 523.0) > 5e-1) {
1452  SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1453  "atom test %d failed: wrong von Mises stress value",
1454  atom_test);
1455  }
1456  if (atom_test == 2 && fabs(von_mises_stress - 16.3) > 5e-2) {
1457  SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1458  "atom test %d failed: wrong von Mises stress value",
1459  atom_test);
1460  }
1461  if (atom_test == 3 && fabs(von_mises_stress - 14.9) > 5e-2) {
1462  SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1463  "atom test %d failed: wrong von Mises stress value",
1464  atom_test);
1465  }
1466  }
1467  }
1468 
1469  MOFEM_LOG_SYNCHRONISE(mField.get_comm());
1470  }
1471 
1472  MoFEMFunctionReturn(0);
1473  };
1474  auto null = boost::shared_ptr<FEMethod>();
1475  CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(), null,
1476  monitor_ptr, null);
1477  MoFEMFunctionReturn(0);
1478  };
1479 
1480  auto set_fieldsplit_preconditioner = [&](auto solver) {
1481  MoFEMFunctionBeginHot;
1482 
1483  SNES snes;
1484  CHKERR TSGetSNES(solver, &snes);
1485  KSP ksp;
1486  CHKERR SNESGetKSP(snes, &ksp);
1487  PC pc;
1488  CHKERR KSPGetPC(ksp, &pc);
1489  PetscBool is_pcfs = PETSC_FALSE;
1490  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1491 
1492  // Setup fieldsplit (block) solver - optional: yes/no
1493  if (is_pcfs == PETSC_TRUE) {
1494  auto bc_mng = mField.getInterface<BcManager>();
1495  auto is_mng = mField.getInterface<ISManager>();
1496  auto name_prb = simple->getProblemName();
1497 
1498  SmartPetscObj<IS> is_u;
1499  CHKERR is_mng->isCreateProblemFieldAndRank(name_prb, ROW, "U", 0,
1500  SPACE_DIM, is_u);
1501  SmartPetscObj<IS> is_flux;
1502  CHKERR is_mng->isCreateProblemFieldAndRank(name_prb, ROW, "FLUX", 0, 0,
1503  is_flux);
1504  SmartPetscObj<IS> is_T;
1505  CHKERR is_mng->isCreateProblemFieldAndRank(name_prb, ROW, "T", 0, 0,
1506  is_T);
1507  SmartPetscObj<IS> is_TBC;
1508  CHKERR is_mng->isCreateProblemFieldAndRank(name_prb, ROW, "TBC", 0, 0,
1509  is_TBC);
1510  IS is_tmp, is_tmp2;
1511  CHKERR ISExpand(is_T, is_flux, &is_tmp);
1512  CHKERR ISExpand(is_TBC, is_tmp, &is_tmp2);
1513  CHKERR ISDestroy(&is_tmp);
1514  auto is_TFlux = SmartPetscObj<IS>(is_tmp2);
1515 
1516  CHKERR ISSort(is_u);
1517  CHKERR ISSort(is_TFlux);
1518  CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_TFlux);
1519  CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_u);
1520  }
1521 
1522  MoFEMFunctionReturnHot(0);
1523  };
1524 
1525  auto D = createDMVector(dm);
1526  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_FORWARD);
1527  CHKERR TSSetSolution(solver, D);
1528  CHKERR TSSetFromOptions(solver);
1529 
1530  CHKERR set_section_monitor(solver);
1531  CHKERR set_fieldsplit_preconditioner(solver);
1532  CHKERR set_time_monitor(dm, solver);
1533 
1534  CHKERR TSSetUp(solver);
1535  CHKERR TSSolve(solver, NULL);
1536 
1537  MoFEMFunctionReturn(0);
1538 }
1539 //! [Solve]
1540 
1541 static char help[] = "...\n\n";
1542 
1543 int main(int argc, char *argv[]) {
1544 
1545  // Initialisation of MoFEM/PETSc and MOAB data structures
1546  const char param_file[] = "param_file.petsc";
1547  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
1548 
1549  // Add logging channel for example
1550  auto core_log = logging::core::get();
1551  core_log->add_sink(
1552  LogManager::createSink(LogManager::getStrmWorld(), "ThermoElastic"));
1553  LogManager::setLog("ThermoElastic");
1554  MOFEM_LOG_TAG("ThermoElastic", "ThermoElastic");
1555 
1556  core_log->add_sink(
1557  LogManager::createSink(LogManager::getStrmSync(), "ThermoElasticSync"));
1558  LogManager::setLog("ThermoElasticSync");
1559  MOFEM_LOG_TAG("ThermoElasticSync", "ThermoElasticSync");
1560 
1561  try {
1562 
1563  //! [Register MoFEM discrete manager in PETSc]
1564  DMType dm_name = "DMMOFEM";
1565  CHKERR DMRegister_MoFEM(dm_name);
1566  //! [Register MoFEM discrete manager in PETSc
1567 
1568  //! [Create MoAB]
1569  moab::Core mb_instance; ///< mesh database
1570  moab::Interface &moab = mb_instance; ///< mesh database interface
1571  //! [Create MoAB]
1572 
1573  //! [Create MoFEM]
1574  MoFEM::Core core(moab); ///< finite element database
1575  MoFEM::Interface &m_field = core; ///< finite element database interface
1576  //! [Create MoFEM]
1577 
1578  //! [Load mesh]
1579  Simple *simple = m_field.getInterface<Simple>();
1580  CHKERR simple->getOptions();
1581  CHKERR simple->loadFile();
1582  //! [Load mesh]
1583 
1584  //! [ThermoElasticProblem]
1585  ThermoElasticProblem ex(m_field);
1586  CHKERR ex.runProblem();
1587  //! [ThermoElasticProblem]
1588  }
1589  CATCH_ERRORS;
1590 
1591  CHKERR MoFEM::Core::Finalize();
1592 }
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
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:589
ThermoElasticProblem::fluxFieldPtr
boost::shared_ptr< MatrixDouble > fluxFieldPtr
Definition: thermo_elastic.cpp:180
ThermoElasticProblem::BlockedParameters::coeffExpansion
VectorDouble coeffExpansion
Definition: thermo_elastic.cpp:196
default_coeff_expansion
double default_coeff_expansion
Definition: thermo_elastic.cpp:132
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
is_plane_strain
PetscBool is_plane_strain
Definition: thermo_elastic.cpp:130
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:1184
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:51
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:168
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
do_output_domain
PetscBool do_output_domain
Definition: thermo_elastic.cpp:146
MoFEM::OpFluxRhsImpl
Definition: Natural.hpp:39
ThermoElasticProblem::BlockedParameters::getHeatCapacityPtr
auto getHeatCapacityPtr()
Definition: thermo_elastic.cpp:213
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
MoFEM::NaturalBC
Natural boundary conditions.
Definition: Natural.hpp:57
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:80
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:209
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
ref_temp
double ref_temp
Definition: thermo_elastic.cpp:127
atom_test
int atom_test
Definition: thermo_elastic.cpp:143
ThermoElasticProblem::fieldEvalCoords
std::array< double, SPACE_DIM > fieldEvalCoords
Definition: thermo_elastic.cpp:176
ThermoElasticProblem::bC
MoFEMErrorCode bC()
[Set up problem]
Definition: thermo_elastic.cpp:641
MoFEM::OpCalculateScalarFieldValuesDot
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
Definition: UserDataOperators.hpp:273
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
ThermalRadiation.hpp
Implementation of thermal radiation bc.
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
ThermoElasticProblem::setupProblem
MoFEMErrorCode setupProblem()
add fields
Definition: thermo_elastic.cpp:558
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
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:126
OpBase
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
FTENSOR_INDEX
#define FTENSOR_INDEX(DIM, I)
Definition: Templates.hpp:2011
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:88
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
sdf.r
int r
Definition: sdf.py:8
ThermoElasticProblem::getCommandLineParameters
MoFEMErrorCode getCommandLineParameters()
[Run problem]
Definition: thermo_elastic.cpp:480
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
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1502
ThermoElasticProblem::tempFieldPtr
boost::shared_ptr< VectorDouble > tempFieldPtr
Definition: thermo_elastic.cpp:179
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
ThermoElasticProblem::dispGradPtr
boost::shared_ptr< MatrixDouble > dispGradPtr
Definition: thermo_elastic.cpp:182
ThermoElasticOps::OpStressThermal
Definition: ThermoElasticOps.hpp:23
ThermoElasticOps.hpp
SPACE_DIM
constexpr int SPACE_DIM
Definition: thermo_elastic.cpp:19
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:133
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:166
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:200
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:94
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
SideEle
ElementsAndOps< SPACE_DIM >::SideEle SideEle
Definition: plastic.cpp:61
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
EntData
EntitiesFieldData::EntData EntData
Definition: thermo_elastic.cpp:22
OpKCauchy
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Linear elastic problem]
Definition: thermo_elastic.cpp:38
ThermoElasticProblem::BlockedParameters::heatCapacity
double heatCapacity
Definition: thermo_elastic.cpp:198
MoFEM::AddFluxToLhsPipelineImpl
Definition: Natural.hpp:49
doEvalField
PetscBool doEvalField
Definition: incompressible_elasticity.cpp:41
ThermalConvection.hpp
Implementation of thermal convection bc.
init_temp
double init_temp
Definition: thermo_elastic.cpp:128
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: ThermalConvection.hpp:12
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
OpInternalForceCauchy
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: thermo_elastic.cpp:42
ThermoElasticProblem::OPs
MoFEMErrorCode OPs()
[Boundary condition]
Definition: thermo_elastic.cpp:766
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', SPACE_DIM >
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:415
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:204
ThermoElasticProblem::strainFieldPtr
boost::shared_ptr< MatrixDouble > strainFieldPtr
Definition: thermo_elastic.cpp:183
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:224
ThermoElasticProblem::dispFieldPtr
boost::shared_ptr< MatrixDouble > dispFieldPtr
Definition: thermo_elastic.cpp:181
order_flux
int order_flux
Definition: thermo_elastic.cpp:140
Range
DomainEleOp
main
int main(int argc, char *argv[])
Definition: thermo_elastic.cpp:1543
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:157
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
shear_modulus_G
double shear_modulus_G
Definition: dynamic_first_order_con_law.cpp:97
save_every
int save_every
Definition: thermo_elastic.cpp:145
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:125
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
MoFEM::EntitiesFieldData::EntData::getFTensor1N
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
Definition: EntitiesFieldData.cpp:640
order_disp
int order_disp
Definition: thermo_elastic.cpp:141
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
ThermoElasticProblem::BlockedParameters
Definition: thermo_elastic.cpp:193
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::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition: PetscSmartObj.hpp:202
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::FormsIntegrators
Integrator forms.
Definition: FormsIntegrators.hpp:306
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
MoFEM::EssentialBC::Assembly
Assembly methods.
Definition: Essential.hpp:119
order_temp
int order_temp
Definition: thermo_elastic.cpp:139
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:67
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:136
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
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:195
ThermoElasticProblem::mField
MoFEM::Interface & mField
Definition: thermo_elastic.cpp:173
DomainEleOp
DomainEle::UserDataOperator DomainEleOp
Definition: thermo_elastic.cpp:24
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:197
ThermoElasticProblem::fieldEvalData
boost::shared_ptr< FieldEvaluatorInterface::SetPtsData > fieldEvalData
Definition: thermo_elastic.cpp:177
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:59
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
OP
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:468
MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition: ForcesAndSourcesCore.hpp:1291
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
do_output_skin
PetscBool do_output_skin
Definition: thermo_elastic.cpp:147
ThermoElasticProblem::stressFieldPtr
boost::shared_ptr< MatrixDouble > stressFieldPtr
Definition: thermo_elastic.cpp:184
help
static char help[]
[Solve]
Definition: thermo_elastic.cpp:1541
ThermoElasticProblem::doEvalField
PetscBool doEvalField
Definition: thermo_elastic.cpp:175
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
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:73
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