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