v0.14.0
seepage.cpp
Go to the documentation of this file.
1 /**
2  * \file seepage.cpp
3  * \example seepage.cpp
4  *
5  * Hydromechanical problem
6  *
7  */
8 
9 /* This file is part of MoFEM.
10  * MoFEM is free software: you can redistribute it and/or modify it under
11  * the terms of the GNU Lesser General Public License as published by the
12  * Free Software Foundation, either version 3 of the License, or (at your
13  * option) any later version.
14  *
15  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
16  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18  * License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public
21  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
22 
23 #ifndef EXECUTABLE_DIMENSION
24 #define EXECUTABLE_DIMENSION 2
25 #endif
26 
27 #include <MoFEM.hpp>
28 #include <MatrixFunction.hpp>
29 #include <IntegrationRules.hpp>
30 
31 
32 using namespace MoFEM;
33 
34 constexpr int SPACE_DIM =
35  EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
36 
37 using EntData = EntitiesFieldData::EntData;
38 using DomainEle = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::DomainEle;
39 using DomainEleOp = DomainEle::UserDataOperator;
40 using BoundaryEle =
41  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
42 using BoundaryEleOp = BoundaryEle::UserDataOperator;
43 using PostProcEle = PostProcBrokenMeshInMoab<DomainEle>;
44 
45 using AssemblyDomainEleOp =
46  FormsIntegrators<DomainEleOp>::Assembly<PETSC>::OpBase;
47 
48 //! [Only used with Hooke equation (linear material model)]
49 using OpKCauchy = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
50  GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
51 using OpInternalForceCauchy = FormsIntegrators<DomainEleOp>::Assembly<
52  PETSC>::LinearForm<GAUSS>::OpGradTimesSymTensor<1, SPACE_DIM, SPACE_DIM>;
53 //! [Only used with Hooke equation (linear material model)]
54 
55 //! [Only used for dynamics]
56 using OpMass = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
57  GAUSS>::OpMass<1, SPACE_DIM>;
58 using OpInertiaForce = FormsIntegrators<DomainEleOp>::Assembly<
59  PETSC>::LinearForm<GAUSS>::OpBaseTimesVector<1, SPACE_DIM, 1>;
60 //! [Only used for dynamics]
61 
62 //! [Only used with Hencky/nonlinear material]
63 using OpKPiola = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
64  GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
65 using OpInternalForcePiola = FormsIntegrators<DomainEleOp>::Assembly<
66  PETSC>::LinearForm<GAUSS>::OpGradTimesTensor<1, SPACE_DIM, SPACE_DIM>;
67 //! [Only used with Hencky/nonlinear material]
68 
69 //! [Essential boundary conditions]
70 using OpBoundaryMass = FormsIntegrators<BoundaryEleOp>::Assembly<
71  PETSC>::BiLinearForm<GAUSS>::OpMass<1, SPACE_DIM>;
72 using OpBoundaryVec = FormsIntegrators<BoundaryEleOp>::Assembly<
73  PETSC>::LinearForm<GAUSS>::OpBaseTimesVector<1, SPACE_DIM, 0>;
74 using OpBoundaryInternal = FormsIntegrators<BoundaryEleOp>::Assembly<
75  PETSC>::LinearForm<GAUSS>::OpBaseTimesVector<1, SPACE_DIM, 1>;
76 //! [Essential boundary conditions]
77 
78 using OpBaseDivU = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
79  GAUSS>::OpMixScalarTimesDiv<SPACE_DIM>;
80 
81 // Thermal operators
82 /**
83  * @brief Integrate Lhs base of flux (1/k) base of flux (FLUX x FLUX)
84  *
85  */
86 using OpHdivHdiv = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
87  GAUSS>::OpMass<3, 3>;
88 
89 /**
90  * @brief Integrate Lhs div of base of flux time base of temperature (FLUX x T)
91  * and transpose of it, i.e. (T x FLAX)
92  *
93  */
94 using OpHdivQ = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
95  GAUSS>::OpMixDivTimesScalar<SPACE_DIM>;
96 
97 /**
98  * @brief Integrate Lhs base of temperature times (heat capacity) times base of
99  * temperature (T x T)
100  *
101  */
102 using OpCapacity = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
103  GAUSS>::OpMass<1, 1>;
104 
105 /**
106  * @brief Integrating Rhs flux base (1/k) flux (FLUX)
107  */
108 using OpHdivFlux = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
109  GAUSS>::OpBaseTimesVector<3, 3, 1>;
110 
111 /**
112  * @brief Integrate Rhs div flux base times temperature (T)
113  *
114  */
115 using OpHDivH = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
116  GAUSS>::OpMixDivTimesU<3, 1, 2>;
117 
118 /**
119  * @brief Integrate Rhs base of temperature time heat capacity times heat rate
120  * (T)
121  *
122  */
123 using OpBaseDotH = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
124  GAUSS>::OpBaseTimesScalarField<1>;
125 
126 /**
127  * @brief Integrate Rhs base of temperature times divergent of flux (T)
128  *
129  */
130 using OpBaseDivFlux = OpBaseDotH;
131 
132 //! [Body and heat source]
133 using DomainNaturalBCRhs =
134  NaturalBC<DomainEleOp>::Assembly<PETSC>::LinearForm<GAUSS>;
135 using OpBodyForce =
136  DomainNaturalBCRhs::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, SPACE_DIM>;
137 using OpHeatSource =
138  DomainNaturalBCRhs::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 1>;
139 using DomainNaturalBCLhs =
140  NaturalBC<DomainEleOp>::Assembly<PETSC>::BiLinearForm<GAUSS>;
141 //! [Body and heat source]
142 
143 //! [Natural boundary conditions]
144 using BoundaryNaturalBC =
145  NaturalBC<BoundaryEleOp>::Assembly<PETSC>::LinearForm<GAUSS>;
146 using OpForce = BoundaryNaturalBC::OpFlux<NaturalForceMeshsets, 1, SPACE_DIM>;
147 using OpTemperatureBC =
148  BoundaryNaturalBC::OpFlux<NaturalTemperatureMeshsets, 3, SPACE_DIM>;
149 //! [Natural boundary conditions]
150 
151 //! [Essential boundary conditions (Least square approach)]
152 using OpEssentialFluxRhs =
153  EssentialBC<BoundaryEleOp>::Assembly<PETSC>::LinearForm<
154  GAUSS>::OpEssentialRhs<HeatFluxCubitBcData, 3, SPACE_DIM>;
155 using OpEssentialFluxLhs =
156  EssentialBC<BoundaryEleOp>::Assembly<PETSC>::BiLinearForm<
157  GAUSS>::OpEssentialLhs<HeatFluxCubitBcData, 3, SPACE_DIM>;
158 //! [Essential boundary conditions (Least square approach)]
159 
160 double scale = 1.;
161 
162 double default_young_modulus = 1;
163 double default_poisson_ratio = 0.25; // zero for verification
164 double default_conductivity = 1;
165 double fluid_density = 10;
166 
167 // #include <OpPostProcElastic.hpp>
168 #include <SeepageOps.hpp>
169 
170 auto save_range = [](moab::Interface &moab, const std::string name,
171  const Range r) {
172  MoFEMFunctionBegin;
173  auto out_meshset = get_temp_meshset_ptr(moab);
174  CHKERR moab.add_entities(*out_meshset, r);
175  CHKERR moab.write_file(name.c_str(), "VTK", "", out_meshset->get_ptr(), 1);
176  MoFEMFunctionReturn(0);
177 };
178 
179 struct Seepage {
180 
181  Seepage(MoFEM::Interface &m_field) : mField(m_field) {}
182 
183  MoFEMErrorCode runProblem();
184 
185 private:
186  MoFEM::Interface &mField;
187 
188  MoFEMErrorCode setupProblem();
189  MoFEMErrorCode createCommonData();
190  MoFEMErrorCode bC();
191  MoFEMErrorCode OPs();
192  MoFEMErrorCode tsSolve();
193 
194  PetscBool doEvalField;
195  std::array<double, SPACE_DIM> fieldEvalCoords;
196 
197  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
198  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
199  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
200 
201  struct BlockedParameters
202  : public boost::enable_shared_from_this<BlockedParameters> {
203  MatrixDouble mD;
204  double fluidConductivity;
205 
206  inline auto getDPtr() {
207  return boost::shared_ptr<MatrixDouble>(shared_from_this(), &mD);
208  }
209 
210  inline auto getConductivityPtr() {
211  return boost::shared_ptr<double>(shared_from_this(), &fluidConductivity);
212  }
213  };
214 
215  MoFEMErrorCode addMatBlockOps(
216  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
217  std::string block_elastic_name, std::string block_thermal_name,
218  boost::shared_ptr<BlockedParameters> blockedParamsPtr, Sev sev);
219 };
220 
221 MoFEMErrorCode Seepage::addMatBlockOps(
222  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
223  std::string block_elastic_name, std::string block_thermal_name,
224  boost::shared_ptr<BlockedParameters> blockedParamsPtr, Sev sev) {
225  MoFEMFunctionBegin;
226 
227  struct OpMatElasticBlocks : public DomainEleOp {
228  OpMatElasticBlocks(boost::shared_ptr<MatrixDouble> m, double bulk_modulus_K,
229  double shear_modulus_G, MoFEM::Interface &m_field,
230  Sev sev,
231  std::vector<const CubitMeshSets *> meshset_vec_ptr)
232  : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), matDPtr(m),
233  bulkModulusKDefault(bulk_modulus_K),
234  shearModulusGDefault(shear_modulus_G) {
235  CHK_THROW_MESSAGE(extractElasticBlockData(m_field, meshset_vec_ptr, sev),
236  "Can not get data from block");
237  }
238 
239  MoFEMErrorCode doWork(int side, EntityType type,
240  EntitiesFieldData::EntData &data) {
241  MoFEMFunctionBegin;
242 
243  for (auto &b : blockData) {
244 
245  if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
246  CHKERR getMatDPtr(matDPtr, b.bulkModulusK, b.shearModulusG);
247  MoFEMFunctionReturnHot(0);
248  }
249  }
250 
251  CHKERR getMatDPtr(matDPtr, bulkModulusKDefault, shearModulusGDefault);
252  MoFEMFunctionReturn(0);
253  }
254 
255  private:
256  boost::shared_ptr<MatrixDouble> matDPtr;
257 
258  struct BlockData {
259  double bulkModulusK;
260  double shearModulusG;
261  Range blockEnts;
262  };
263 
264  double bulkModulusKDefault;
265  double shearModulusGDefault;
266  std::vector<BlockData> blockData;
267 
268  MoFEMErrorCode
269  extractElasticBlockData(MoFEM::Interface &m_field,
270  std::vector<const CubitMeshSets *> meshset_vec_ptr,
271  Sev sev) {
272  MoFEMFunctionBegin;
273 
274  for (auto m : meshset_vec_ptr) {
275  MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Elastic Block") << *m;
276  std::vector<double> block_data;
277  CHKERR m->getAttributes(block_data);
278  if (block_data.size() < 2) {
279  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
280  "Expected that block has two attributes");
281  }
282  auto get_block_ents = [&]() {
283  Range ents;
284  CHKERR
285  m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
286  return ents;
287  };
288 
289  double young_modulus = block_data[0];
290  double poisson_ratio = block_data[1];
291  double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
292  double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
293 
294  MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Elastic Block")
295  << m->getName() << ": E = " << young_modulus
296  << " nu = " << poisson_ratio;
297 
298  blockData.push_back(
299  {bulk_modulus_K, shear_modulus_G, get_block_ents()});
300  }
301  MOFEM_LOG_CHANNEL("WORLD");
302  MoFEMFunctionReturn(0);
303  }
304 
305  MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
306  double bulk_modulus_K, double shear_modulus_G) {
307  MoFEMFunctionBegin;
308  //! [Calculate elasticity tensor]
309  auto set_material_stiffness = [&]() {
310  FTensor::Index<'i', SPACE_DIM> i;
311  FTensor::Index<'j', SPACE_DIM> j;
312  FTensor::Index<'k', SPACE_DIM> k;
313  FTensor::Index<'l', SPACE_DIM> l;
314  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
315  double A = (SPACE_DIM == 2)
316  ? 2 * shear_modulus_G /
317  (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
318  : 1;
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_bulk_modulus_K, 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 OpMatFluidBlocks : public DomainEleOp {
352  OpMatFluidBlocks(boost::shared_ptr<double> conductivity_ptr,
353  MoFEM::Interface &m_field, Sev sev,
354  std::vector<const CubitMeshSets *> meshset_vec_ptr)
355  : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE),
356  conductivityPtr(conductivity_ptr) {
357  CHK_THROW_MESSAGE(extractThermalBlockData(m_field, meshset_vec_ptr, sev),
358  "Can not get data from block");
359  }
360 
361  MoFEMErrorCode doWork(int side, EntityType type,
362  EntitiesFieldData::EntData &data) {
363  MoFEMFunctionBegin;
364 
365  for (auto &b : blockData) {
366 
367  if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
368  *conductivityPtr = b.conductivity;
369  MoFEMFunctionReturnHot(0);
370  }
371  }
372 
373  *conductivityPtr = default_conductivity;
374 
375  MoFEMFunctionReturn(0);
376  }
377 
378  private:
379  struct BlockData {
380  double conductivity;
381  Range blockEnts;
382  };
383 
384  std::vector<BlockData> blockData;
385 
386  MoFEMErrorCode
387  extractThermalBlockData(MoFEM::Interface &m_field,
388  std::vector<const CubitMeshSets *> meshset_vec_ptr,
389  Sev sev) {
390  MoFEMFunctionBegin;
391 
392  for (auto m : meshset_vec_ptr) {
393  MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Thermal Block") << *m;
394  std::vector<double> block_data;
395  CHKERR m->getAttributes(block_data);
396  if (block_data.size() < 1) {
397  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
398  "Expected that block has two attributes");
399  }
400  auto get_block_ents = [&]() {
401  Range ents;
402  CHKERR
403  m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
404  return ents;
405  };
406 
407  MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Thermal Block")
408  << m->getName() << ": conductivity = " << block_data[0];
409 
410  blockData.push_back({block_data[0], get_block_ents()});
411  }
412  MOFEM_LOG_CHANNEL("WORLD");
413  MoFEMFunctionReturn(0);
414  }
415 
416  boost::shared_ptr<double> expansionPtr;
417  boost::shared_ptr<double> conductivityPtr;
418  boost::shared_ptr<double> capacityPtr;
419  };
420 
421  pipeline.push_back(new OpMatFluidBlocks(
422  blockedParamsPtr->getConductivityPtr(), mField, sev,
423 
424  // Get blockset using regular expression
425  mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
426 
427  (boost::format("%s(.*)") % block_thermal_name).str()
428 
429  ))
430 
431  ));
432 
433  MoFEMFunctionReturn(0);
434 }
435 
436 //! [Run problem]
437 MoFEMErrorCode Seepage::runProblem() {
438  MoFEMFunctionBegin;
439  CHKERR setupProblem();
440  CHKERR createCommonData();
441  CHKERR bC();
442  CHKERR OPs();
443  CHKERR tsSolve();
444  MoFEMFunctionReturn(0);
445 }
446 //! [Run problem]
447 
448 //! [Set up problem]
449 MoFEMErrorCode Seepage::setupProblem() {
450  MoFEMFunctionBegin;
451  Simple *simple = mField.getInterface<Simple>();
452  // Add field
453  constexpr FieldApproximationBase base = AINSWORTH_LEGENDRE_BASE;
454  // Mechanical fields
455  CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
456  CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
457  // Temerature
458  const auto flux_space = (SPACE_DIM == 2) ? HCURL : HDIV;
459  CHKERR simple->addDomainField("H", L2, AINSWORTH_LEGENDRE_BASE, 1);
460  CHKERR simple->addDomainField("FLUX", flux_space, DEMKOWICZ_JACOBI_BASE, 1);
461  CHKERR simple->addBoundaryField("FLUX", flux_space, DEMKOWICZ_JACOBI_BASE, 1);
462 
463  int order = 2.;
464  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
465  CHKERR simple->setFieldOrder("U", order);
466  CHKERR simple->setFieldOrder("H", order - 1);
467  CHKERR simple->setFieldOrder("FLUX", order);
468 
469  CHKERR simple->setUp();
470 
471  MoFEMFunctionReturn(0);
472 }
473 //! [Set up problem]
474 
475 //! [Create common data]
476 MoFEMErrorCode Seepage::createCommonData() {
477  MoFEMFunctionBegin;
478 
479  auto get_command_line_parameters = [&]() {
480  MoFEMFunctionBegin;
481  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
482  &default_young_modulus, PETSC_NULL);
483  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
484  &default_poisson_ratio, PETSC_NULL);
485  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-conductivity",
486  &default_conductivity, PETSC_NULL);
487  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-fluid_density",
488  &fluid_density, PETSC_NULL);
489 
490  MOFEM_LOG("Seepage", Sev::inform)
491  << "Default Young modulus " << default_young_modulus;
492  MOFEM_LOG("Seepage", Sev::inform)
493  << "Default Poisson ratio " << default_poisson_ratio;
494  MOFEM_LOG("Seepage", Sev::inform)
495  << "Default conductivity " << default_conductivity;
496  MOFEM_LOG("Seepage", Sev::inform) << "Fluid denisty " << fluid_density;
497 
498  int coords_dim = SPACE_DIM;
499  CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
500  fieldEvalCoords.data(), &coords_dim,
501  &doEvalField);
502 
503  MoFEMFunctionReturn(0);
504  };
505 
506  CHKERR get_command_line_parameters();
507 
508  MoFEMFunctionReturn(0);
509 }
510 //! [Create common data]
511 
512 //! [Boundary condition]
513 MoFEMErrorCode Seepage::bC() {
514  MoFEMFunctionBegin;
515 
516  MOFEM_LOG("SYNC", Sev::noisy) << "bC";
517  MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::noisy);
518 
519  auto simple = mField.getInterface<Simple>();
520  auto bc_mng = mField.getInterface<BcManager>();
521 
522  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
523  simple->getProblemName(), "U");
524  CHKERR bc_mng->pushMarkDOFsOnEntities<HeatFluxCubitBcData>(
525  simple->getProblemName(), "FLUX", false);
526 
527  auto get_skin = [&]() {
528  Range body_ents;
529  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
530  Skinner skin(&mField.get_moab());
531  Range skin_ents;
532  CHKERR skin.find_skin(0, body_ents, false, skin_ents);
533  return skin_ents;
534  };
535 
536  auto filter_flux_blocks = [&](auto skin) {
537  auto remove_cubit_blocks = [&](auto c) {
538  MoFEMFunctionBegin;
539  for (auto m :
540 
541  mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(c)
542 
543  ) {
544  Range ents;
545  CHKERR mField.get_moab().get_entities_by_dimension(
546  m->getMeshset(), SPACE_DIM - 1, ents, true);
547  skin = subtract(skin, ents);
548  }
549  MoFEMFunctionReturn(0);
550  };
551 
552  auto remove_named_blocks = [&](auto n) {
553  MoFEMFunctionBegin;
554  for (auto m : mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
555  std::regex(
556 
557  (boost::format("%s(.*)") % n).str()
558 
559  ))
560 
561  ) {
562  Range ents;
563  CHKERR mField.get_moab().get_entities_by_dimension(
564  m->getMeshset(), SPACE_DIM - 1, ents, true);
565  skin = subtract(skin, ents);
566  }
567  MoFEMFunctionReturn(0);
568  };
569 
570  CHK_THROW_MESSAGE(remove_cubit_blocks(NODESET | TEMPERATURESET),
571  "remove_cubit_blocks");
572  CHK_THROW_MESSAGE(remove_cubit_blocks(SIDESET | HEATFLUXSET),
573  "remove_cubit_blocks");
574  CHK_THROW_MESSAGE(remove_named_blocks("PRESSURE"), "remove_named_blocks");
575  CHK_THROW_MESSAGE(remove_named_blocks("FLUIDFLUX"), "remove_named_blocks");
576 
577  return skin;
578  };
579 
580  auto filter_true_skin = [&](auto skin) {
581  Range boundary_ents;
582  ParallelComm *pcomm =
583  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
584  CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
585  PSTATUS_NOT, -1, &boundary_ents);
586  return boundary_ents;
587  };
588 
589  auto remove_flux_ents = filter_true_skin(filter_flux_blocks(get_skin()));
590 
591  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
592  remove_flux_ents);
593 
594  MOFEM_LOG("SYNC", Sev::noisy) << remove_flux_ents << endl;
595  MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::noisy);
596 
597 #ifndef NDEBUG
598 
599  CHKERR save_range(
600  mField.get_moab(),
601  (boost::format("flux_remove_%d.vtk") % mField.get_comm_rank()).str(),
602  remove_flux_ents);
603 
604 #endif
605 
606  CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
607  simple->getProblemName(), "FLUX", remove_flux_ents);
608 
609  MoFEMFunctionReturn(0);
610 }
611 //! [Boundary condition]
612 
613 //! [Push operators to pipeline]
614 MoFEMErrorCode Seepage::OPs() {
615  MoFEMFunctionBegin;
616  auto pipeline_mng = mField.getInterface<PipelineManager>();
617  auto simple = mField.getInterface<Simple>();
618  auto bc_mng = mField.getInterface<BcManager>();
619 
620  auto boundary_marker =
621  bc_mng->getMergedBlocksMarker(vector<string>{"FLUIDFLUX"});
622 
623  auto u_grad_ptr = boost::make_shared<MatrixDouble>();
624  auto dot_u_grad_ptr = boost::make_shared<MatrixDouble>();
625  auto trace_dot_u_grad_ptr = boost::make_shared<VectorDouble>();
626  auto h_ptr = boost::make_shared<VectorDouble>();
627  auto dot_h_ptr = boost::make_shared<VectorDouble>();
628  auto flux_ptr = boost::make_shared<MatrixDouble>();
629  auto div_flux_ptr = boost::make_shared<VectorDouble>();
630  auto strain_ptr = boost::make_shared<MatrixDouble>();
631  auto stress_ptr = boost::make_shared<MatrixDouble>();
632 
633  auto time_scale = boost::make_shared<TimeScale>();
634 
635  auto block_params = boost::make_shared<BlockedParameters>();
636  auto mDPtr = block_params->getDPtr();
637  auto conductivity_ptr = block_params->getConductivityPtr();
638 
639  auto integration_rule = [](int, int, int approx_order) {
640  return 2 * approx_order;
641  };
642 
643  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
644  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
645  CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
646  CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
647 
648  auto add_domain_base_ops = [&](auto &pip) {
649  MoFEMFunctionBegin;
650 
651  CHKERR addMatBlockOps(pip, "MAT_ELASTIC", "MAT_FLUID", block_params,
652  Sev::inform);
653  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV});
654 
655  pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
656  "U", u_grad_ptr));
657  pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(u_grad_ptr, strain_ptr));
658  pip.push_back(new OpCalculateVectorFieldGradientDot<SPACE_DIM, SPACE_DIM>(
659  "U", dot_u_grad_ptr));
660  pip.push_back(new OpCalculateTraceFromMat<SPACE_DIM>(dot_u_grad_ptr,
661  trace_dot_u_grad_ptr));
662 
663  pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
664  pip.push_back(new OpCalculateScalarFieldValuesDot("H", dot_h_ptr));
665  pip.push_back(new OpCalculateHVecVectorField<3>("FLUX", flux_ptr));
666  pip.push_back(new OpCalculateHdivVectorDivergence<3, SPACE_DIM>(
667  "FLUX", div_flux_ptr));
668 
669  MoFEMFunctionReturn(0);
670  };
671 
672  auto add_domain_ops_lhs_mechanical = [&](auto &pip) {
673  MoFEMFunctionBegin;
674  pip.push_back(new OpKCauchy("U", "U", mDPtr));
675  pip.push_back(new OpBaseDivU(
676  "H", "U",
677  [](const double, const double, const double) { return -9.81; }, true,
678  true));
679  MoFEMFunctionReturn(0);
680  };
681 
682  auto add_domain_ops_rhs_mechanical = [&](auto &pip) {
683  MoFEMFunctionBegin;
684 
685  CHKERR DomainNaturalBCRhs::AddFluxToPipeline<OpBodyForce>::add(
686  pip, mField, "U", {time_scale}, "BODY_FORCE", Sev::inform);
687 
688  // Calculate internal forece
689  pip.push_back(new OpTensorTimesSymmetricTensor<SPACE_DIM, SPACE_DIM>(
690  "U", strain_ptr, stress_ptr, mDPtr));
691  pip.push_back(new OpInternalForceCauchy("U", stress_ptr));
692  pip.push_back(
693  new SeepageOps::OpDomainRhsHydrostaticStress<SPACE_DIM>("U", h_ptr));
694 
695  MoFEMFunctionReturn(0);
696  };
697 
698  auto add_domain_ops_lhs_seepage = [&](auto &pip, auto &fe) {
699  MoFEMFunctionBegin;
700  auto resistance = [conductivity_ptr](const double, const double,
701  const double) {
702  return (1. / *(conductivity_ptr));
703  };
704 
705  auto unity = []() constexpr { return -1; };
706  pip.push_back(new OpHdivHdiv("FLUX", "FLUX", resistance));
707  pip.push_back(new OpHdivQ("FLUX", "H", unity, true));
708  auto op_base_div_u = new OpBaseDivU(
709  "H", "U", [](double, double, double) constexpr { return -1; }, false,
710  false);
711  op_base_div_u->feScalingFun = [](const FEMethod *fe_ptr) {
712  return fe_ptr->ts_a;
713  };
714  pip.push_back(op_base_div_u);
715 
716  MoFEMFunctionReturn(0);
717  };
718 
719  auto add_domain_ops_rhs_seepage = [&](auto &pip) {
720  MoFEMFunctionBegin;
721  auto resistance = [conductivity_ptr](double, double, double) {
722  return (1. / (*conductivity_ptr));
723  };
724  auto minus_one = [](double, double, double) constexpr { return -1; };
725 
726  pip.push_back(new OpHdivFlux("FLUX", flux_ptr, resistance));
727  pip.push_back(new OpHDivH("FLUX", h_ptr, minus_one));
728  pip.push_back(new OpBaseDotH("H", trace_dot_u_grad_ptr, minus_one));
729  pip.push_back(new OpBaseDivFlux("H", div_flux_ptr, minus_one));
730 
731  MoFEMFunctionReturn(0);
732  };
733 
734  auto add_boundary_rhs_ops = [&](auto &pip) {
735  MoFEMFunctionBegin;
736 
737  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {HDIV});
738 
739  pip.push_back(new OpSetBc("FLUX", true, boundary_marker));
740 
741  CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpForce>::add(
742  pipeline_mng->getOpBoundaryRhsPipeline(), mField, "U", {time_scale},
743  "FORCE", Sev::inform);
744 
745  CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpTemperatureBC>::add(
746  pipeline_mng->getOpBoundaryRhsPipeline(), mField, "FLUX", {time_scale},
747  "PRESSURE", Sev::inform);
748 
749  pip.push_back(new OpUnSetBc("FLUX"));
750 
751  auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
752  pip.push_back(
753  new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
754  CHKERR EssentialBC<BoundaryEleOp>::Assembly<PETSC>::LinearForm<GAUSS>::
755  AddEssentialToPipeline<OpEssentialFluxRhs>::add(
756  mField, pip, simple->getProblemName(), "FLUX", mat_flux_ptr,
757  {time_scale});
758 
759  MoFEMFunctionReturn(0);
760  };
761 
762  auto add_boundary_lhs_ops = [&](auto &pip) {
763  MoFEMFunctionBegin;
764 
765  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {HDIV});
766 
767  CHKERR EssentialBC<BoundaryEleOp>::Assembly<PETSC>::BiLinearForm<GAUSS>::
768  AddEssentialToPipeline<OpEssentialFluxLhs>::add(
769  mField, pip, simple->getProblemName(), "FLUX");
770 
771  MoFEMFunctionReturn(0);
772  };
773 
774  // LHS
775  CHKERR add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
776  CHKERR add_domain_ops_lhs_mechanical(pipeline_mng->getOpDomainLhsPipeline());
777  CHKERR add_domain_ops_lhs_seepage(pipeline_mng->getOpDomainLhsPipeline(),
778  pipeline_mng->getDomainLhsFE());
779 
780  // RHS
781  CHKERR add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
782  CHKERR add_domain_ops_rhs_mechanical(pipeline_mng->getOpDomainRhsPipeline());
783  CHKERR add_domain_ops_rhs_seepage(pipeline_mng->getOpDomainRhsPipeline());
784 
785  CHKERR add_boundary_rhs_ops(pipeline_mng->getOpBoundaryRhsPipeline());
786  CHKERR add_boundary_lhs_ops(pipeline_mng->getOpBoundaryLhsPipeline());
787 
788  MoFEMFunctionReturn(0);
789 }
790 //! [Push operators to pipeline]
791 
792 //! [Solve]
793 MoFEMErrorCode Seepage::tsSolve() {
794  MoFEMFunctionBegin;
795  auto simple = mField.getInterface<Simple>();
796  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
797  ISManager *is_manager = mField.getInterface<ISManager>();
798 
799  auto dm = simple->getDM();
800  auto solver = pipeline_mng->createTSIM();
801  auto snes_ctx_ptr = getDMSnesCtx(dm);
802 
803  auto set_section_monitor = [&](auto solver) {
804  MoFEMFunctionBegin;
805  SNES snes;
806  CHKERR TSGetSNES(solver, &snes);
807  CHKERR SNESMonitorSet(snes,
808  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
809  void *))MoFEMSNESMonitorFields,
810  (void *)(snes_ctx_ptr.get()), nullptr);
811  MoFEMFunctionReturn(0);
812  };
813 
814  auto create_post_process_element = [&]() {
815  auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
816 
817  auto block_params = boost::make_shared<BlockedParameters>();
818  auto mDPtr = block_params->getDPtr();
819  CHKERR addMatBlockOps(post_proc_fe->getOpPtrVector(), "MAT_ELASTIC",
820  "MAT_FLUID", block_params, Sev::verbose);
821  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
822  post_proc_fe->getOpPtrVector(), {H1, HDIV});
823 
824  auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
825  auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
826  auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
827 
828  auto h_ptr = boost::make_shared<VectorDouble>();
829  auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
830 
831  post_proc_fe->getOpPtrVector().push_back(
832  new OpCalculateScalarFieldValues("H", h_ptr));
833  post_proc_fe->getOpPtrVector().push_back(
834  new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
835 
836  auto u_ptr = boost::make_shared<MatrixDouble>();
837  post_proc_fe->getOpPtrVector().push_back(
838  new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
839  post_proc_fe->getOpPtrVector().push_back(
840  new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>("U",
841  mat_grad_ptr));
842  post_proc_fe->getOpPtrVector().push_back(
843  new OpSymmetrizeTensor<SPACE_DIM>(mat_grad_ptr, mat_strain_ptr));
844  post_proc_fe->getOpPtrVector().push_back(
845  new OpTensorTimesSymmetricTensor<SPACE_DIM, SPACE_DIM>(
846  "U", mat_strain_ptr, mat_stress_ptr, mDPtr));
847 
848  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
849 
850  post_proc_fe->getOpPtrVector().push_back(
851 
852  new OpPPMap(
853 
854  post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
855 
856  {{"H", h_ptr}},
857 
858  {{"U", u_ptr}, {"FLUX", mat_flux_ptr}},
859 
860  {},
861 
862  {{"STRAIN", mat_strain_ptr}, {"STRESS", mat_stress_ptr}}
863 
864  )
865 
866  );
867 
868  return post_proc_fe;
869  };
870 
871  auto create_creaction_fe = [&]() {
872  auto fe_ptr = boost::make_shared<DomainEle>(mField);
873  fe_ptr->getRuleHook = [](int, int, int o) { return 2 * o; };
874 
875  auto &pip = fe_ptr->getOpPtrVector();
876 
877  auto block_params = boost::make_shared<BlockedParameters>();
878  auto mDPtr = block_params->getDPtr();
879  CHKERR addMatBlockOps(pip, "MAT_ELASTIC", "MAT_FLUID", block_params,
880  Sev::verbose);
881  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV});
882 
883  auto u_grad_ptr = boost::make_shared<MatrixDouble>();
884  auto strain_ptr = boost::make_shared<MatrixDouble>();
885  auto stress_ptr = boost::make_shared<MatrixDouble>();
886 
887  pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
888  "U", u_grad_ptr));
889  pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(u_grad_ptr, strain_ptr));
890 
891  // Calculate internal forece
892  pip.push_back(new OpTensorTimesSymmetricTensor<SPACE_DIM, SPACE_DIM>(
893  "U", strain_ptr, stress_ptr, mDPtr));
894  pip.push_back(new OpInternalForceCauchy("U", stress_ptr));
895 
896  fe_ptr->postProcessHook =
897  EssentialPreProcReaction<DisplacementCubitBcData>(mField, fe_ptr);
898 
899  return fe_ptr;
900  };
901 
902  auto monitor_ptr = boost::make_shared<FEMethod>();
903  auto post_proc_fe = create_post_process_element();
904  auto res = createDMVector(dm);
905  auto rections_fe = create_creaction_fe();
906 
907  auto set_time_monitor = [&](auto dm, auto solver) {
908  MoFEMFunctionBegin;
909  monitor_ptr->preProcessHook = [&]() {
910  MoFEMFunctionBegin;
911 
912  CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
913  post_proc_fe,
914  monitor_ptr->getCacheWeakPtr());
915  CHKERR post_proc_fe->writeFile(
916  "out_" + boost::lexical_cast<std::string>(monitor_ptr->ts_step) +
917  ".h5m");
918 
919  rections_fe->f = res;
920  CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
921  rections_fe,
922  monitor_ptr->getCacheWeakPtr());
923 
924  double nrm;
925  CHKERR VecNorm(res, NORM_2, &nrm);
926  MOFEM_LOG("Seepage", Sev::verbose)
927  << "Residual norm " << nrm << " at time step "
928  << monitor_ptr->ts_step;
929 
930  if (doEvalField) {
931 
932  auto scalar_field_ptr = boost::make_shared<VectorDouble>();
933  auto vector_field_ptr = boost::make_shared<MatrixDouble>();
934  auto tensor_field_ptr = boost::make_shared<MatrixDouble>();
935 
936  auto field_eval_data = mField.getInterface<FieldEvaluatorInterface>()
937  ->getData<DomainEle>();
938 
939  if constexpr (SPACE_DIM == 3) {
940  CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree3D(
941  field_eval_data, simple->getDomainFEName());
942  } else {
943  CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree2D(
944  field_eval_data, simple->getDomainFEName());
945  }
946 
947  field_eval_data->setEvalPoints(fieldEvalCoords.data(), 1);
948  auto no_rule = [](int, int, int) { return -1; };
949 
950  auto field_eval_ptr = field_eval_data->feMethodPtr.lock();
951  field_eval_ptr->getRuleHook = no_rule;
952 
953  auto u_grad_ptr = boost::make_shared<MatrixDouble>();
954  auto strain_ptr = boost::make_shared<MatrixDouble>();
955  auto stress_ptr = boost::make_shared<MatrixDouble>();
956  auto h_ptr = boost::make_shared<VectorDouble>();
957 
958  auto block_params = boost::make_shared<BlockedParameters>();
959  auto mDPtr = block_params->getDPtr();
960  CHKERR addMatBlockOps(field_eval_ptr->getOpPtrVector(), "MAT_ELASTIC",
961  "MAT_FLUID", block_params, Sev::noisy);
962  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
963  field_eval_ptr->getOpPtrVector(), {H1, HDIV});
964  field_eval_ptr->getOpPtrVector().push_back(
965  new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
966  "U", u_grad_ptr));
967  field_eval_ptr->getOpPtrVector().push_back(
968  new OpSymmetrizeTensor<SPACE_DIM>(u_grad_ptr, strain_ptr));
969  field_eval_ptr->getOpPtrVector().push_back(
970  new OpCalculateScalarFieldValues("H", h_ptr));
971  field_eval_ptr->getOpPtrVector().push_back(
972  new OpTensorTimesSymmetricTensor<SPACE_DIM, SPACE_DIM>(
973  "U", strain_ptr, stress_ptr, mDPtr));
974 
975  if constexpr (SPACE_DIM == 3) {
976  CHKERR mField.getInterface<FieldEvaluatorInterface>()
977  ->evalFEAtThePoint3D(
978  fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
979  simple->getDomainFEName(), field_eval_data,
980  mField.get_comm_rank(), mField.get_comm_rank(), nullptr,
981  MF_EXIST, QUIET);
982  } else {
983  CHKERR mField.getInterface<FieldEvaluatorInterface>()
984  ->evalFEAtThePoint2D(
985  fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
986  simple->getDomainFEName(), field_eval_data,
987  mField.get_comm_rank(), mField.get_comm_rank(), nullptr,
988  MF_EXIST, QUIET);
989  }
990 
991  MOFEM_LOG("SeepageSync", Sev::inform)
992  << "Eval point hydrostatic hight: " << *h_ptr;
993  MOFEM_LOG("SeepageSync", Sev::inform)
994  << "Eval point skeleton stress pressure: " << *stress_ptr;
995  MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::inform);
996  }
997 
998  MoFEMFunctionReturn(0);
999  };
1000  auto null = boost::shared_ptr<FEMethod>();
1001  CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(), null,
1002  monitor_ptr, null);
1003  MoFEMFunctionReturn(0);
1004  };
1005 
1006  auto set_fieldsplit_preconditioner = [&](auto solver) {
1007  MoFEMFunctionBeginHot;
1008 
1009  SNES snes;
1010  CHKERR TSGetSNES(solver, &snes);
1011  KSP ksp;
1012  CHKERR SNESGetKSP(snes, &ksp);
1013  PC pc;
1014  CHKERR KSPGetPC(ksp, &pc);
1015  PetscBool is_pcfs = PETSC_FALSE;
1016  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1017 
1018  // Setup fieldsplit (block) solver - optional: yes/no
1019  if (is_pcfs == PETSC_TRUE) {
1020  auto bc_mng = mField.getInterface<BcManager>();
1021  auto is_mng = mField.getInterface<ISManager>();
1022  auto name_prb = simple->getProblemName();
1023 
1024  SmartPetscObj<IS> is_u;
1025  CHKERR is_mng->isCreateProblemFieldAndRank(name_prb, ROW, "U", 0,
1026  SPACE_DIM, is_u);
1027  SmartPetscObj<IS> is_flux;
1028  CHKERR is_mng->isCreateProblemFieldAndRank(name_prb, ROW, "FLUX", 0, 0,
1029  is_flux);
1030  SmartPetscObj<IS> is_H;
1031  CHKERR is_mng->isCreateProblemFieldAndRank(name_prb, ROW, "H", 0, 0,
1032  is_H);
1033  IS is_tmp;
1034  CHKERR ISExpand(is_H, is_flux, &is_tmp);
1035  auto is_Flux = SmartPetscObj<IS>(is_tmp);
1036 
1037  auto is_all_bc = bc_mng->getBlockIS(name_prb, "FLUIDFLUX", "FLUX", 0, 0);
1038  int is_all_bc_size;
1039  CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
1040  MOFEM_LOG("ThermoElastic", Sev::inform)
1041  << "Field split block size " << is_all_bc_size;
1042  if (is_all_bc_size) {
1043  IS is_tmp2;
1044  CHKERR ISDifference(is_Flux, is_all_bc, &is_tmp2);
1045  is_Flux = SmartPetscObj<IS>(is_tmp2);
1046  CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
1047  is_all_bc); // boundary block
1048  }
1049 
1050  CHKERR ISSort(is_u);
1051  CHKERR ISSort(is_Flux);
1052  CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_Flux);
1053  CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_u);
1054  }
1055 
1056  MoFEMFunctionReturnHot(0);
1057  };
1058 
1059  auto pre_proc_ptr = boost::make_shared<FEMethod>();
1060  auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
1061  auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
1062  auto time_scale = boost::make_shared<TimeScale>();
1063 
1064  auto get_bc_hook_rhs = [this, pre_proc_ptr, time_scale]() {
1065  EssentialPreProc<DisplacementCubitBcData> hook(mField, pre_proc_ptr,
1066  {time_scale}, false);
1067  return hook;
1068  };
1069 
1070  auto get_post_proc_hook_rhs = [this, post_proc_rhs_ptr]() {
1071  MoFEMFunctionBegin;
1072  CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
1073  mField, post_proc_rhs_ptr, nullptr, Sev::verbose)();
1074  CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(
1075  mField, post_proc_rhs_ptr, 1.)();
1076  MoFEMFunctionReturn(0);
1077  };
1078  auto get_post_proc_hook_lhs = [this, post_proc_lhs_ptr]() {
1079  return EssentialPostProcLhs<DisplacementCubitBcData>(mField,
1080  post_proc_lhs_ptr, 1.);
1081  };
1082 
1083  pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
1084  post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
1085  post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
1086 
1087  auto ts_ctx_ptr = getDMTsCtx(dm);
1088  ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
1089  ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
1090  ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
1091  ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
1092 
1093  auto D = createDMVector(dm);
1094  CHKERR TSSetSolution(solver, D);
1095  CHKERR TSSetFromOptions(solver);
1096 
1097  CHKERR set_section_monitor(solver);
1098  CHKERR set_fieldsplit_preconditioner(solver);
1099  CHKERR set_time_monitor(dm, solver);
1100 
1101  CHKERR TSSetUp(solver);
1102  CHKERR TSSolve(solver, NULL);
1103 
1104  MoFEMFunctionReturn(0);
1105 }
1106 //! [Solve]
1107 
1108 static char help[] = "...\n\n";
1109 
1110 int main(int argc, char *argv[]) {
1111 
1112  // Initialisation of MoFEM/PETSc and MOAB data structures
1113  const char param_file[] = "param_file.petsc";
1114  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
1115 
1116  // Add logging channel for example
1117  auto core_log = logging::core::get();
1118  core_log->add_sink(
1119  LogManager::createSink(LogManager::getStrmWorld(), "Seepage"));
1120  LogManager::setLog("Seepage");
1121  MOFEM_LOG_TAG("Seepage", "seepage");
1122  core_log->add_sink(
1123  LogManager::createSink(LogManager::getStrmSync(), "SeepageSync"));
1124  LogManager::setLog("SeepageSync");
1125  MOFEM_LOG_TAG("SeepageSync", "seepage");
1126 
1127  try {
1128 
1129  //! [Register MoFEM discrete manager in PETSc]
1130  DMType dm_name = "DMMOFEM";
1131  CHKERR DMRegister_MoFEM(dm_name);
1132  //! [Register MoFEM discrete manager in PETSc
1133 
1134  //! [Create MoAB]
1135  moab::Core mb_instance; ///< mesh database
1136  moab::Interface &moab = mb_instance; ///< mesh database interface
1137  //! [Create MoAB]
1138 
1139  //! [Create MoFEM]
1140  MoFEM::Core core(moab); ///< finite element database
1141  MoFEM::Interface &m_field = core; ///< finite element database insterface
1142  //! [Create MoFEM]
1143 
1144  //! [Load mesh]
1145  Simple *simple = m_field.getInterface<Simple>();
1146  CHKERR simple->getOptions();
1147  CHKERR simple->loadFile();
1148  //! [Load mesh]
1149 
1150  //! [Seepage]
1151  Seepage ex(m_field);
1152  CHKERR ex.runProblem();
1153  //! [Seepage]
1154  }
1155  CATCH_ERRORS;
1156 
1157  CHKERR MoFEM::Core::Finalize();
1158 }
Seepage::BlockedParameters::fluidConductivity
double fluidConductivity
Definition: seepage.cpp:204
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
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
SIDESET
@ SIDESET
Definition: definitions.h:160
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
save_range
auto save_range
Definition: seepage.cpp:170
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
MoFEM::EssentialPreProcReaction< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:157
Seepage::createCommonData
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: seepage.cpp:476
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
SeepageOps.hpp
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
help
static char help[]
[Solve]
Definition: seepage.cpp:1108
young_modulus
double young_modulus
Young modulus.
Definition: plastic.cpp:121
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
default_conductivity
double default_conductivity
Definition: seepage.cpp:164
SPACE_DIM
constexpr int SPACE_DIM
Definition: seepage.cpp:34
MoFEM::PipelineManager::ElementsAndOpsByDim
Definition: PipelineManager.hpp:38
EntData
EntitiesFieldData::EntData EntData
Definition: seepage.cpp:37
MoFEM::EssentialPostProcLhs< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:134
DomainEleOp
DomainEle::UserDataOperator DomainEleOp
Definition: seepage.cpp:39
MoFEM::NaturalBC::Assembly::BiLinearForm
Definition: Natural.hpp:74
Seepage::OPs
MoFEMErrorCode OPs()
[Boundary condition]
Definition: seepage.cpp:614
OpKCauchy
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Only used with Hooke equation (linear material model)]
Definition: seepage.cpp:50
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
Seepage::doEvalField
PetscBool doEvalField
Definition: seepage.cpp:194
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
OpHdivFlux
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 3, 3, 1 > OpHdivFlux
Integrating Rhs flux base (1/k) flux (FLUX)
Definition: seepage.cpp:109
MoFEM::EssentialBC
Essential boundary conditions.
Definition: Essential.hpp:111
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
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::OpCalculateScalarFieldValuesDot
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
Definition: UserDataOperators.hpp:273
Seepage::uZScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: seepage.cpp:199
MoFEM::EssentialPostProcRhs< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:113
OpKPiola
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: seepage.cpp:64
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
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
Seepage::mField
MoFEM::Interface & mField
Definition: seepage.cpp:186
MoFEM::DisplacementCubitBcData
Definition of the displacement bc data structure.
Definition: BCData.hpp:76
OpHdivHdiv
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, 3 > OpHdivHdiv
Integrate Lhs base of flux (1/k) base of flux (FLUX x FLUX)
Definition: seepage.cpp:87
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Seepage::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: seepage.cpp:221
OpBaseDotH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpBaseDotH
Integrate Rhs base of temperature time heat capacity times heat rate (T)
Definition: seepage.cpp:124
OpBase
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
OpInternalForceCauchy
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: seepage.cpp:52
sdf.r
int r
Definition: sdf.py:8
Seepage::BlockedParameters::mD
MatrixDouble mD
Definition: seepage.cpp:203
Seepage::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition: seepage.cpp:449
Seepage::bC
MoFEMErrorCode bC()
[Create common data]
Definition: seepage.cpp:513
order
constexpr int order
Definition: dg_projection.cpp:18
OpMixScalarTimesDiv
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMixScalarTimesDiv< SPACE_DIM, COORD_TYPE > OpMixScalarTimesDiv
Definition: free_surface.cpp:125
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
OpBoundaryMass
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: seepage.cpp:71
ROW
@ ROW
Definition: definitions.h:136
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
scale
double scale
[Essential boundary conditions (Least square approach)]
Definition: seepage.cpp:160
OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
Definition: dynamic_first_order_con_law.cpp:63
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:667
NODESET
@ NODESET
Definition: definitions.h:159
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:178
OpHDivH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 1, 2 > OpHDivH
Integrate Rhs div flux base times temperature (T)
Definition: seepage.cpp:116
default_young_modulus
double default_young_modulus
Definition: seepage.cpp:162
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
OpInternalForcePiola
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: seepage.cpp:66
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
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
Seepage::BlockedParameters::getConductivityPtr
auto getConductivityPtr()
Definition: seepage.cpp:210
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
double
MoFEM::OpEssentialLhsImpl
Enforce essential constrains on lhs.
Definition: Essential.hpp:81
convert.type
type
Definition: convert.py:64
Seepage
Definition: seepage.cpp:179
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:317
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
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::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
EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition: seepage.cpp:24
MatrixFunction.hpp
OpKCauchy
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Linear elastic problem]
Definition: thermo_elastic.cpp:34
doEvalField
PetscBool doEvalField
Definition: incompressible_elasticity.cpp:41
Seepage::uXScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: seepage.cpp:197
Seepage::tsSolve
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
Definition: seepage.cpp:793
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
size_symm
constexpr auto size_symm
Definition: plastic.cpp:42
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::OpTensorTimesSymmetricTensor
Calculate .
Definition: UserDataOperators.hpp:1883
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
BiLinearForm
OpBoundaryInternal
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
Definition: seepage.cpp:75
default_poisson_ratio
double default_poisson_ratio
Definition: seepage.cpp:163
OpInternalForceCauchy
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: thermo_elastic.cpp:38
OpGradTimesTensor
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
Definition: operators_tests.cpp:48
Seepage::BlockedParameters::getDPtr
auto getDPtr()
Definition: seepage.cpp:206
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
Seepage::fieldEvalCoords
std::array< double, SPACE_DIM > fieldEvalCoords
Definition: seepage.cpp:195
fluid_density
double fluid_density
Definition: seepage.cpp:165
Seepage::Seepage
Seepage(MoFEM::Interface &m_field)
Definition: seepage.cpp:181
integration_rule
auto integration_rule
Definition: free_surface.cpp:185
Seepage::runProblem
MoFEMErrorCode runProblem()
[Run problem]
Definition: seepage.cpp:437
MoFEM::OpUnSetBc
Definition: FormsIntegrators.hpp:49
Range
DomainEleOp
OpBaseDivFlux
OpBaseDotH OpBaseDivFlux
Integrate Rhs base of temperature times divergent of flux (T)
Definition: seepage.cpp:130
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
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
OpHdivQ
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< SPACE_DIM > OpHdivQ
Integrate Lhs div of base of flux time base of temperature (FLUX x T) and transpose of it,...
Definition: seepage.cpp:95
shear_modulus_G
double shear_modulus_G
Definition: dynamic_first_order_con_law.cpp:97
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::OpSymmetrizeTensor
Definition: UserDataOperators.hpp:1948
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
IntegrationRules.hpp
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
Seepage::BlockedParameters
Definition: seepage.cpp:201
approx_order
int approx_order
Definition: test_broken_space.cpp:50
MoFEM::OpCalculateTraceFromMat
Calculates the trace of an input matrix.
Definition: UserDataOperators.hpp:3319
MoFEM::DMMoFEMTSSetMonitor
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition: DMMoFEM.cpp:1056
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::OpSetBc
Set indices on entities on finite element.
Definition: FormsIntegrators.hpp:38
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
TEMPERATURESET
@ TEMPERATURESET
Definition: definitions.h:168
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::BlockData
Definition: MeshsetsManager.cpp:755
MoFEM::OpCalculateHdivVectorDivergence
Calculate divergence of vector field.
Definition: UserDataOperators.hpp:2212
MoFEM::MoFEMSNESMonitorFields
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
Definition: SnesCtx.cpp:232
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
MoFEM::EssentialBC::Assembly
Assembly methods.
Definition: Essential.hpp:119
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
main
int main(int argc, char *argv[])
Definition: seepage.cpp:1110
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
MoFEM::getDMSnesCtx
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition: DMMoFEM.hpp:1127
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
HEATFLUXSET
@ HEATFLUXSET
Definition: definitions.h:169
QUIET
@ QUIET
Definition: definitions.h:221
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
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::SmartPetscObj< VecScatter >
OpBoundaryVec
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
Definition: seepage.cpp:73
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
MF_EXIST
@ MF_EXIST
Definition: definitions.h:113
convert.int
int
Definition: convert.py:64
MoFEM::OpCalculateVectorFieldGradientDot
Get field gradients time derivative at integration pts for scalar filed rank 0, i....
Definition: UserDataOperators.hpp:1550
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::getDMTsCtx
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1141
Seepage::uYScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: seepage.cpp:198
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
OpBaseDivU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixScalarTimesDiv< SPACE_DIM > OpBaseDivU
[Essential boundary conditions]
Definition: seepage.cpp:79
SeepageOps::OpDomainRhsHydrostaticStress
Definition: SeepageOps.hpp:22
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