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