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