v0.15.0
Loading...
Searching...
No Matches
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
31using namespace MoFEM;
32
33constexpr int SPACE_DIM =
34 EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
35
38using DomainEleOp = DomainEle::UserDataOperator;
41using BoundaryEleOp = BoundaryEle::UserDataOperator;
43
46
47//! [Only used with Hooke equation (linear material model)]
49 GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
52//! [Only used with Hooke equation (linear material model)]
53
54//! [Only used for dynamics]
59//! [Only used for dynamics]
60
61//! [Only used with Hencky/nonlinear material]
63 GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
66//! [Only used with Hencky/nonlinear material]
67
68//! [Essential boundary conditions]
75//! [Essential boundary conditions]
76
79
80// Thermal operators
81/**
82 * @brief Integrate Lhs base of flux (1/k) base of flux (FLUX x FLUX)
83 *
84 */
87
88/**
89 * @brief Integrate Lhs div of base of flux time base of temperature (FLUX x T)
90 * and transpose of it, i.e. (T x FLAX)
91 *
92 */
94 GAUSS>::OpMixDivTimesScalar<SPACE_DIM>;
95
96/**
97 * @brief Integrate Lhs base of temperature times (heat capacity) times base of
98 * temperature (T x T)
99 *
100 */
103
104/**
105 * @brief Integrating Rhs flux base (1/k) flux (FLUX)
106 */
108 GAUSS>::OpBaseTimesVector<3, 3, 1>;
109
110/**
111 * @brief Integrate Rhs div flux base times temperature (T)
112 *
113 */
115 GAUSS>::OpMixDivTimesU<3, 1, 2>;
116
117/**
118 * @brief Integrate Rhs base of temperature time heat capacity times heat rate
119 * (T)
120 *
121 */
123 GAUSS>::OpBaseTimesScalarField<1>;
124
125/**
126 * @brief Integrate Rhs base of temperature times divergent of flux (T)
127 *
128 */
130
131//! [Body and heat source]
140//! [Body and heat source]
141
142//! [Natural boundary conditions]
148//! [Natural boundary conditions]
149
150//! [Essential boundary conditions (Least square approach)]
153 GAUSS>::OpEssentialRhs<HeatFluxCubitBcData, 3, SPACE_DIM>;
156 GAUSS>::OpEssentialLhs<HeatFluxCubitBcData, 3, SPACE_DIM>;
157//! [Essential boundary conditions (Least square approach)]
158
159double scale = 1.;
160
162double default_poisson_ratio = 0.25; // zero for verification
164double fluid_density = 10;
165
166// #include <OpPostProcElastic.hpp>
167#include <SeepageOps.hpp>
168
169auto save_range = [](moab::Interface &moab, const std::string name,
170 const Range r) {
172 auto out_meshset = get_temp_meshset_ptr(moab);
173 CHKERR moab.add_entities(*out_meshset, r);
174 CHKERR moab.write_file(name.c_str(), "VTK", "", out_meshset->get_ptr(), 1);
176};
177
178struct Seepage {
179
180 Seepage(MoFEM::Interface &m_field) : mField(m_field) {}
181
183
184private:
186
192
193 PetscBool doEvalField;
194 std::array<double, SPACE_DIM> fieldEvalCoords;
195
196 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
197 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
198 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
199
201 : public boost::enable_shared_from_this<BlockedParameters> {
204
205 inline auto getDPtr() {
206 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &mD);
207 }
208
209 inline auto getConductivityPtr() {
210 return boost::shared_ptr<double>(shared_from_this(), &fluidConductivity);
211 }
212 };
213
215 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
216 std::string block_elastic_name, std::string block_thermal_name,
217 boost::shared_ptr<BlockedParameters> blockedParamsPtr, Sev sev);
218};
219
221 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
222 std::string block_elastic_name, std::string block_thermal_name,
223 boost::shared_ptr<BlockedParameters> blockedParamsPtr, Sev sev) {
225
226 struct OpMatElasticBlocks : public DomainEleOp {
227 OpMatElasticBlocks(boost::shared_ptr<MatrixDouble> m, double bulk_modulus_K,
228 double shear_modulus_G, MoFEM::Interface &m_field,
229 Sev sev,
230 std::vector<const CubitMeshSets *> meshset_vec_ptr)
231 : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), matDPtr(m),
232 bulkModulusKDefault(bulk_modulus_K),
233 shearModulusGDefault(shear_modulus_G) {
234 CHK_THROW_MESSAGE(extractElasticBlockData(m_field, meshset_vec_ptr, sev),
235 "Can not get data from block");
236 }
237
238 MoFEMErrorCode doWork(int side, EntityType type,
241
242 for (auto &b : blockData) {
243
244 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
245 CHKERR getMatDPtr(matDPtr, b.bulkModulusK, b.shearModulusG);
247 }
248 }
249
250 CHKERR getMatDPtr(matDPtr, bulkModulusKDefault, shearModulusGDefault);
252 }
253
254 private:
255 boost::shared_ptr<MatrixDouble> matDPtr;
256
257 struct BlockData {
258 double bulkModulusK;
259 double shearModulusG;
260 Range blockEnts;
261 };
262
263 double bulkModulusKDefault;
264 double shearModulusGDefault;
265 std::vector<BlockData> blockData;
266
268 extractElasticBlockData(MoFEM::Interface &m_field,
269 std::vector<const CubitMeshSets *> meshset_vec_ptr,
270 Sev sev) {
272
273 for (auto m : meshset_vec_ptr) {
274 MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Elastic Block") << *m;
275 std::vector<double> block_data;
276 CHKERR m->getAttributes(block_data);
277 if (block_data.size() < 2) {
278 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
279 "Expected that block has two attributes");
280 }
281 auto get_block_ents = [&]() {
282 Range ents;
283 CHKERR
284 m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
285 return ents;
286 };
287
288 double young_modulus = block_data[0];
289 double poisson_ratio = block_data[1];
290 double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
291 double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
292
293 MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Elastic Block")
294 << m->getName() << ": E = " << young_modulus
295 << " nu = " << poisson_ratio;
296
297 blockData.push_back(
298 {bulk_modulus_K, shear_modulus_G, get_block_ents()});
299 }
300 MOFEM_LOG_CHANNEL("WORLD");
302 }
303
304 MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
305 double bulk_modulus_K, double shear_modulus_G) {
307 //! [Calculate elasticity tensor]
308 auto set_material_stiffness = [&]() {
314 double A = (SPACE_DIM == 2)
315 ? 2 * shear_modulus_G /
316 (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
317 : 1;
319 t_D(i, j, k, l) =
320 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
321 A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) * t_kd(i, j) *
322 t_kd(k, l);
323 };
324 //! [Calculate elasticity tensor]
325 constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
326 mat_D_ptr->resize(size_symm * size_symm, 1);
327 set_material_stiffness();
329 }
330 };
331
332 double default_bulk_modulus_K =
334 double default_shear_modulus_G =
336
337 pipeline.push_back(new OpMatElasticBlocks(
338 blockedParamsPtr->getDPtr(), default_bulk_modulus_K,
339 default_shear_modulus_G, mField, sev,
340
341 // Get blockset using regular expression
343
344 (boost::format("%s(.*)") % block_elastic_name).str()
345
346 ))
347
348 ));
349
350 struct OpMatFluidBlocks : public DomainEleOp {
351 OpMatFluidBlocks(boost::shared_ptr<double> conductivity_ptr,
352 MoFEM::Interface &m_field, Sev sev,
353 std::vector<const CubitMeshSets *> meshset_vec_ptr)
354 : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE),
355 conductivityPtr(conductivity_ptr) {
356 CHK_THROW_MESSAGE(extractThermalBlockData(m_field, meshset_vec_ptr, sev),
357 "Can not get data from block");
358 }
359
360 MoFEMErrorCode doWork(int side, EntityType type,
363
364 for (auto &b : blockData) {
365
366 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
367 *conductivityPtr = b.conductivity;
369 }
370 }
371
372 *conductivityPtr = default_conductivity;
373
375 }
376
377 private:
378 struct BlockData {
379 double conductivity;
380 Range blockEnts;
381 };
382
383 std::vector<BlockData> blockData;
384
386 extractThermalBlockData(MoFEM::Interface &m_field,
387 std::vector<const CubitMeshSets *> meshset_vec_ptr,
388 Sev sev) {
390
391 for (auto m : meshset_vec_ptr) {
392 MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Thermal Block") << *m;
393 std::vector<double> block_data;
394 CHKERR m->getAttributes(block_data);
395 if (block_data.size() < 1) {
396 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
397 "Expected that block has two attributes");
398 }
399 auto get_block_ents = [&]() {
400 Range ents;
401 CHKERR
402 m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
403 return ents;
404 };
405
406 MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Thermal Block")
407 << m->getName() << ": conductivity = " << block_data[0];
408
409 blockData.push_back({block_data[0], get_block_ents()});
410 }
411 MOFEM_LOG_CHANNEL("WORLD");
413 }
414
415 boost::shared_ptr<double> expansionPtr;
416 boost::shared_ptr<double> conductivityPtr;
417 boost::shared_ptr<double> capacityPtr;
418 };
419
420 pipeline.push_back(new OpMatFluidBlocks(
421 blockedParamsPtr->getConductivityPtr(), mField, sev,
422
423 // Get blockset using regular expression
425
426 (boost::format("%s(.*)") % block_thermal_name).str()
427
428 ))
429
430 ));
431
433}
434
435//! [Run problem]
445//! [Run problem]
446
447//! [Set up problem]
451 // Add field
453 // Mechanical fields
455 CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
456 // Temerature
457 const auto flux_space = (SPACE_DIM == 2) ? HCURL : HDIV;
458 CHKERR simple->addDomainField("H", L2, AINSWORTH_LEGENDRE_BASE, 1);
459 CHKERR simple->addDomainField("FLUX", flux_space, DEMKOWICZ_JACOBI_BASE, 1);
460 CHKERR simple->addBoundaryField("FLUX", flux_space, DEMKOWICZ_JACOBI_BASE, 1);
461
462 int order = 2.;
463 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order, PETSC_NULLPTR);
464 CHKERR simple->setFieldOrder("U", order);
465 CHKERR simple->setFieldOrder("H", order - 1);
466 CHKERR simple->setFieldOrder("FLUX", order);
467
468 CHKERR simple->setUp();
469
471}
472//! [Set up problem]
473
474//! [Create common data]
477
478 auto get_command_line_parameters = [&]() {
480 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-young_modulus",
481 &default_young_modulus, PETSC_NULLPTR);
482 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-poisson_ratio",
483 &default_poisson_ratio, PETSC_NULLPTR);
484 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-conductivity",
485 &default_conductivity, PETSC_NULLPTR);
486 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-fluid_density",
487 &fluid_density, PETSC_NULLPTR);
488
489 MOFEM_LOG("Seepage", Sev::inform)
490 << "Default Young modulus " << default_young_modulus;
491 MOFEM_LOG("Seepage", Sev::inform)
492 << "Default Poisson ratio " << default_poisson_ratio;
493 MOFEM_LOG("Seepage", Sev::inform)
494 << "Default conductivity " << default_conductivity;
495 MOFEM_LOG("Seepage", Sev::inform) << "Fluid denisty " << fluid_density;
496
497 int coords_dim = 3;
498 CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
499 fieldEvalCoords.data(), &coords_dim,
500 &doEvalField);
501
503 };
504
505 CHKERR get_command_line_parameters();
506
508}
509//! [Create common data]
510
511//! [Boundary condition]
514
515 MOFEM_LOG("SYNC", Sev::noisy) << "bC";
517
519 auto bc_mng = mField.getInterface<BcManager>();
520
521 CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
522 simple->getProblemName(), "U");
523 CHKERR bc_mng->pushMarkDOFsOnEntities<HeatFluxCubitBcData>(
524 simple->getProblemName(), "FLUX", false);
525
526 auto get_skin = [&]() {
527 Range body_ents;
528 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
529 Skinner skin(&mField.get_moab());
530 Range skin_ents;
531 CHKERR skin.find_skin(0, body_ents, false, skin_ents);
532 return skin_ents;
533 };
534
535 auto filter_flux_blocks = [&](auto skin) {
536 auto remove_cubit_blocks = [&](auto c) {
538 for (auto m :
539
541
542 ) {
543 Range ents;
544 CHKERR mField.get_moab().get_entities_by_dimension(
545 m->getMeshset(), SPACE_DIM - 1, ents, true);
546 skin = subtract(skin, ents);
547 }
549 };
550
551 auto remove_named_blocks = [&](auto n) {
554 std::regex(
555
556 (boost::format("%s(.*)") % n).str()
557
558 ))
559
560 ) {
561 Range ents;
562 CHKERR mField.get_moab().get_entities_by_dimension(
563 m->getMeshset(), SPACE_DIM - 1, ents, true);
564 skin = subtract(skin, ents);
565 }
567 };
568
569 CHK_THROW_MESSAGE(remove_cubit_blocks(NODESET | TEMPERATURESET),
570 "remove_cubit_blocks");
571 CHK_THROW_MESSAGE(remove_cubit_blocks(SIDESET | HEATFLUXSET),
572 "remove_cubit_blocks");
573 CHK_THROW_MESSAGE(remove_named_blocks("PRESSURE"), "remove_named_blocks");
574 CHK_THROW_MESSAGE(remove_named_blocks("FLUIDFLUX"), "remove_named_blocks");
575
576 return skin;
577 };
578
579 auto filter_true_skin = [&](auto skin) {
580 Range boundary_ents;
581 ParallelComm *pcomm =
582 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
583 CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
584 PSTATUS_NOT, -1, &boundary_ents);
585 return boundary_ents;
586 };
587
588 auto remove_flux_ents = filter_true_skin(filter_flux_blocks(get_skin()));
589
590 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
591 remove_flux_ents);
592
593 MOFEM_LOG("SYNC", Sev::noisy) << remove_flux_ents << endl;
595
596#ifndef NDEBUG
597
600 (boost::format("flux_remove_%d.vtk") % mField.get_comm_rank()).str(),
601 remove_flux_ents);
602
603#endif
604
605 CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
606 simple->getProblemName(), "FLUX", remove_flux_ents);
607
609}
610//! [Boundary condition]
611
612//! [Push operators to pipeline]
615 auto pipeline_mng = mField.getInterface<PipelineManager>();
617 auto bc_mng = mField.getInterface<BcManager>();
618
619 auto boundary_marker =
620 bc_mng->getMergedBlocksMarker(vector<string>{"FLUIDFLUX"});
621
622 auto u_grad_ptr = boost::make_shared<MatrixDouble>();
623 auto dot_u_grad_ptr = boost::make_shared<MatrixDouble>();
624 auto trace_dot_u_grad_ptr = boost::make_shared<VectorDouble>();
625 auto h_ptr = boost::make_shared<VectorDouble>();
626 auto dot_h_ptr = boost::make_shared<VectorDouble>();
627 auto flux_ptr = boost::make_shared<MatrixDouble>();
628 auto div_flux_ptr = boost::make_shared<VectorDouble>();
629 auto strain_ptr = boost::make_shared<MatrixDouble>();
630 auto stress_ptr = boost::make_shared<MatrixDouble>();
631
632 auto time_scale = boost::make_shared<TimeScale>();
633
634 auto block_params = boost::make_shared<BlockedParameters>();
635 auto mDPtr = block_params->getDPtr();
636 auto conductivity_ptr = block_params->getConductivityPtr();
637
638 auto integration_rule = [](int, int, int approx_order) {
639 return 2 * approx_order;
640 };
641
642 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
643 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
644 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
645 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
646
647 auto add_domain_base_ops = [&](auto &pip) {
649
650 CHKERR addMatBlockOps(pip, "MAT_ELASTIC", "MAT_FLUID", block_params,
651 Sev::inform);
653
655 "U", u_grad_ptr));
656 pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(u_grad_ptr, strain_ptr));
658 "U", dot_u_grad_ptr));
659 pip.push_back(new OpCalculateTraceFromMat<SPACE_DIM>(dot_u_grad_ptr,
660 trace_dot_u_grad_ptr));
661
662 pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
663 pip.push_back(new OpCalculateScalarFieldValuesDot("H", dot_h_ptr));
664 pip.push_back(new OpCalculateHVecVectorField<3>("FLUX", flux_ptr));
666 "FLUX", div_flux_ptr));
667
669 };
670
671 auto add_domain_ops_lhs_mechanical = [&](auto &pip) {
673 pip.push_back(new OpKCauchy("U", "U", mDPtr));
674 pip.push_back(new OpBaseDivU(
675 "H", "U",
676 [](const double, const double, const double) { return -9.81; }, true,
677 true));
679 };
680
681 auto add_domain_ops_rhs_mechanical = [&](auto &pip) {
683
685 pip, mField, "U", {time_scale}, "BODY_FORCE", Sev::inform);
686
687 // Calculate internal forece
689 strain_ptr, stress_ptr, mDPtr));
690 pip.push_back(new OpInternalForceCauchy("U", stress_ptr));
691 pip.push_back(
693
695 };
696
697 auto add_domain_ops_lhs_seepage = [&](auto &pip, auto &fe) {
699 auto resistance = [conductivity_ptr](const double, const double,
700 const double) {
701 return (1. / *(conductivity_ptr));
702 };
703
704 auto unity = []() constexpr { return -1; };
705 pip.push_back(new OpHdivHdiv("FLUX", "FLUX", resistance));
706 pip.push_back(new OpHdivQ("FLUX", "H", unity, true));
707 auto op_base_div_u = new OpBaseDivU(
708 "H", "U", [](double, double, double) constexpr { return -1; }, false,
709 false);
710 op_base_div_u->feScalingFun = [](const FEMethod *fe_ptr) {
711 return fe_ptr->ts_a;
712 };
713 pip.push_back(op_base_div_u);
714
716 };
717
718 auto add_domain_ops_rhs_seepage = [&](auto &pip) {
720 auto resistance = [conductivity_ptr](double, double, double) {
721 return (1. / (*conductivity_ptr));
722 };
723 auto minus_one = [](double, double, double) constexpr { return -1; };
724
725 pip.push_back(new OpHdivFlux("FLUX", flux_ptr, resistance));
726 pip.push_back(new OpHDivH("FLUX", h_ptr, minus_one));
727 pip.push_back(new OpBaseDotH("H", trace_dot_u_grad_ptr, minus_one));
728 pip.push_back(new OpBaseDivFlux("H", div_flux_ptr, minus_one));
729
731 };
732
733 auto add_boundary_rhs_ops = [&](auto &pip) {
735
737
738 pip.push_back(new OpSetBc("FLUX", true, boundary_marker));
739
741 pipeline_mng->getOpBoundaryRhsPipeline(), mField, "U", {time_scale},
742 "FORCE", Sev::inform);
743
745 pipeline_mng->getOpBoundaryRhsPipeline(), mField, "FLUX", {time_scale},
746 "PRESSURE", Sev::inform);
747
748 pip.push_back(new OpUnSetBc("FLUX"));
749
750 auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
751 pip.push_back(
752 new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
755 mField, pip, simple->getProblemName(), "FLUX", mat_flux_ptr,
756 {time_scale});
757
759 };
760
761 auto add_boundary_lhs_ops = [&](auto &pip) {
763
765
768 mField, pip, simple->getProblemName(), "FLUX");
769
771 };
772
773 // LHS
774 CHKERR add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
775 CHKERR add_domain_ops_lhs_mechanical(pipeline_mng->getOpDomainLhsPipeline());
776 CHKERR add_domain_ops_lhs_seepage(pipeline_mng->getOpDomainLhsPipeline(),
777 pipeline_mng->getDomainLhsFE());
778
779 // RHS
780 CHKERR add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
781 CHKERR add_domain_ops_rhs_mechanical(pipeline_mng->getOpDomainRhsPipeline());
782 CHKERR add_domain_ops_rhs_seepage(pipeline_mng->getOpDomainRhsPipeline());
783
784 CHKERR add_boundary_rhs_ops(pipeline_mng->getOpBoundaryRhsPipeline());
785 CHKERR add_boundary_lhs_ops(pipeline_mng->getOpBoundaryLhsPipeline());
786
788}
789//! [Push operators to pipeline]
790
791//! [Solve]
796
797 auto dm = simple->getDM();
798 auto solver = pipeline_mng->createTSIM();
799 auto snes_ctx_ptr = getDMSnesCtx(dm);
800
801 auto set_section_monitor = [&](auto solver) {
803 SNES snes;
804 CHKERR TSGetSNES(solver, &snes);
805 CHKERR SNESMonitorSet(snes,
806 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
808 (void *)(snes_ctx_ptr.get()), nullptr);
810 };
811
812 auto create_post_process_element = [&]() {
813 auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
814
815 auto block_params = boost::make_shared<BlockedParameters>();
816 auto mDPtr = block_params->getDPtr();
817 CHKERR addMatBlockOps(post_proc_fe->getOpPtrVector(), "MAT_ELASTIC",
818 "MAT_FLUID", block_params, Sev::verbose);
820 post_proc_fe->getOpPtrVector(), {H1, HDIV});
821
822 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
823 auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
824 auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
825
826 auto h_ptr = boost::make_shared<VectorDouble>();
827 auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
828
829 post_proc_fe->getOpPtrVector().push_back(
830 new OpCalculateScalarFieldValues("H", h_ptr));
831 post_proc_fe->getOpPtrVector().push_back(
832 new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
833
834 auto u_ptr = boost::make_shared<MatrixDouble>();
835 post_proc_fe->getOpPtrVector().push_back(
837 post_proc_fe->getOpPtrVector().push_back(
839 mat_grad_ptr));
840 post_proc_fe->getOpPtrVector().push_back(
841 new OpSymmetrizeTensor<SPACE_DIM>(mat_grad_ptr, mat_strain_ptr));
842 post_proc_fe->getOpPtrVector().push_back(
844 mat_strain_ptr, mat_stress_ptr, mDPtr));
845
847
848 post_proc_fe->getOpPtrVector().push_back(
849
850 new OpPPMap(
851
852 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
853
854 {{"H", h_ptr}},
855
856 {{"U", u_ptr}, {"FLUX", mat_flux_ptr}},
857
858 {},
859
860 {{"STRAIN", mat_strain_ptr}, {"STRESS", mat_stress_ptr}}
861
862 )
863
864 );
865
866 return post_proc_fe;
867 };
868
869 auto create_creaction_fe = [&]() {
870 auto fe_ptr = boost::make_shared<DomainEle>(mField);
871 fe_ptr->getRuleHook = [](int, int, int o) { return 2 * o; };
872
873 auto &pip = fe_ptr->getOpPtrVector();
874
875 auto block_params = boost::make_shared<BlockedParameters>();
876 auto mDPtr = block_params->getDPtr();
877 CHKERR addMatBlockOps(pip, "MAT_ELASTIC", "MAT_FLUID", block_params,
878 Sev::verbose);
880
881 auto u_grad_ptr = boost::make_shared<MatrixDouble>();
882 auto strain_ptr = boost::make_shared<MatrixDouble>();
883 auto stress_ptr = boost::make_shared<MatrixDouble>();
884
886 "U", u_grad_ptr));
887 pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(u_grad_ptr, strain_ptr));
888
889 // Calculate internal forece
891 strain_ptr, stress_ptr, mDPtr));
892 pip.push_back(new OpInternalForceCauchy("U", stress_ptr));
893
894 fe_ptr->postProcessHook =
896
897 return fe_ptr;
898 };
899
900 auto monitor_ptr = boost::make_shared<FEMethod>();
901 auto post_proc_fe = create_post_process_element();
902 auto res = createDMVector(dm);
903 auto rections_fe = create_creaction_fe();
904
905 auto set_time_monitor = [&](auto dm, auto solver) {
907 monitor_ptr->preProcessHook = [&]() {
909
910 CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
911 post_proc_fe,
912 monitor_ptr->getCacheWeakPtr());
913 CHKERR post_proc_fe->writeFile(
914 "out_" + boost::lexical_cast<std::string>(monitor_ptr->ts_step) +
915 ".h5m");
916
917 rections_fe->f = res;
918 CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
919 rections_fe,
920 monitor_ptr->getCacheWeakPtr());
921
922 CHKERR VecAssemblyBegin(res);
923 CHKERR VecAssemblyEnd(res);
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 CHKERR mField.getInterface<FieldEvaluatorInterface>()
940 ->buildTree<SPACE_DIM>(field_eval_data, simple->getDomainFEName());
941
942 field_eval_data->setEvalPoints(fieldEvalCoords.data(), 1);
943 auto no_rule = [](int, int, int) { return -1; };
944
945 auto field_eval_ptr = field_eval_data->feMethodPtr.lock();
946 field_eval_ptr->getRuleHook = no_rule;
947
948 auto u_grad_ptr = boost::make_shared<MatrixDouble>();
949 auto strain_ptr = boost::make_shared<MatrixDouble>();
950 auto stress_ptr = boost::make_shared<MatrixDouble>();
951 auto h_ptr = boost::make_shared<VectorDouble>();
952
953 auto block_params = boost::make_shared<BlockedParameters>();
954 auto mDPtr = block_params->getDPtr();
955 CHKERR addMatBlockOps(field_eval_ptr->getOpPtrVector(), "MAT_ELASTIC",
956 "MAT_FLUID", block_params, Sev::noisy);
958 field_eval_ptr->getOpPtrVector(), {H1, HDIV});
959 field_eval_ptr->getOpPtrVector().push_back(
961 "U", u_grad_ptr));
962 field_eval_ptr->getOpPtrVector().push_back(
963 new OpSymmetrizeTensor<SPACE_DIM>(u_grad_ptr, strain_ptr));
964 field_eval_ptr->getOpPtrVector().push_back(
965 new OpCalculateScalarFieldValues("H", h_ptr));
966 field_eval_ptr->getOpPtrVector().push_back(
968 strain_ptr, stress_ptr, mDPtr));
969
970 CHKERR mField.getInterface<FieldEvaluatorInterface>()
971 ->evalFEAtThePoint<SPACE_DIM>(
972 fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
973 simple->getDomainFEName(), field_eval_data,
974 mField.get_comm_rank(), mField.get_comm_rank(), nullptr,
975 MF_EXIST, QUIET);
976
977 MOFEM_LOG("SeepageSync", Sev::inform)
978 << "Eval point hydrostatic hight: " << *h_ptr;
979 MOFEM_LOG("SeepageSync", Sev::inform)
980 << "Eval point skeleton stress pressure: " << *stress_ptr;
981 MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::inform);
982 }
983
985 };
986 auto null = boost::shared_ptr<FEMethod>();
987 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(), null,
988 monitor_ptr, null);
990 };
991
992 auto set_fieldsplit_preconditioner = [&](auto solver) {
994
995 SNES snes;
996 CHKERR TSGetSNES(solver, &snes);
997 KSP ksp;
998 CHKERR SNESGetKSP(snes, &ksp);
999 PC pc;
1000 CHKERR KSPGetPC(ksp, &pc);
1001 PetscBool is_pcfs = PETSC_FALSE;
1002 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1003
1004 // Setup fieldsplit (block) solver - optional: yes/no
1005 if (is_pcfs == PETSC_TRUE) {
1006 auto bc_mng = mField.getInterface<BcManager>();
1007 auto is_mng = mField.getInterface<ISManager>();
1008 auto name_prb = simple->getProblemName();
1009
1010 SmartPetscObj<IS> is_u;
1011 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb, ROW, "U", 0,
1012 SPACE_DIM, is_u);
1013 SmartPetscObj<IS> is_flux;
1014 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb, ROW, "FLUX", 0, 0,
1015 is_flux);
1016 SmartPetscObj<IS> is_H;
1017 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb, ROW, "H", 0, 0,
1018 is_H);
1019 IS is_tmp;
1020 CHKERR ISExpand(is_H, is_flux, &is_tmp);
1021 auto is_Flux = SmartPetscObj<IS>(is_tmp);
1022
1023 auto is_all_bc = bc_mng->getBlockIS(name_prb, "FLUIDFLUX", "FLUX", 0, 0);
1024 int is_all_bc_size;
1025 CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
1026 MOFEM_LOG("ThermoElastic", Sev::inform)
1027 << "Field split block size " << is_all_bc_size;
1028 if (is_all_bc_size) {
1029 IS is_tmp2;
1030 CHKERR ISDifference(is_Flux, is_all_bc, &is_tmp2);
1031 is_Flux = SmartPetscObj<IS>(is_tmp2);
1032 CHKERR PCFieldSplitSetIS(pc, PETSC_NULLPTR,
1033 is_all_bc); // boundary block
1034 }
1035
1036 CHKERR ISSort(is_u);
1037 CHKERR ISSort(is_Flux);
1038 CHKERR PCFieldSplitSetIS(pc, PETSC_NULLPTR, is_Flux);
1039 CHKERR PCFieldSplitSetIS(pc, PETSC_NULLPTR, is_u);
1040 }
1041
1043 };
1044
1045 auto pre_proc_ptr = boost::make_shared<FEMethod>();
1046 auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
1047 auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
1048 auto time_scale = boost::make_shared<TimeScale>();
1049
1050 auto get_bc_hook_rhs = [this, pre_proc_ptr, time_scale]() {
1051 EssentialPreProc<DisplacementCubitBcData> hook(mField, pre_proc_ptr,
1052 {time_scale}, false);
1053 return hook;
1054 };
1055
1056 auto get_post_proc_hook_rhs = [this, post_proc_rhs_ptr]() {
1059 mField, post_proc_rhs_ptr, nullptr, Sev::verbose)();
1061 mField, post_proc_rhs_ptr, 1.)();
1063 };
1064 auto get_post_proc_hook_lhs = [this, post_proc_lhs_ptr]() {
1066 post_proc_lhs_ptr, 1.);
1067 };
1068
1069 pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
1070 post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
1071 post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
1072
1073 auto ts_ctx_ptr = getDMTsCtx(dm);
1074 ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
1075 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
1076 ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
1077 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
1078
1079 auto B = createDMMatrix(dm);
1080 CHKERR TSSetIJacobian(solver, B, B, PETSC_NULLPTR, PETSC_NULLPTR);
1081 auto D = createDMVector(dm);
1082 CHKERR TSSetSolution(solver, D);
1083 CHKERR TSSetFromOptions(solver);
1084
1085 CHKERR set_section_monitor(solver);
1086 CHKERR set_fieldsplit_preconditioner(solver);
1087 CHKERR set_time_monitor(dm, solver);
1088
1089 CHKERR TSSetUp(solver);
1090 CHKERR TSSolve(solver, NULL);
1091
1093}
1094//! [Solve]
1095
1096static char help[] = "...\n\n";
1097
1098int main(int argc, char *argv[]) {
1099
1100 // Initialisation of MoFEM/PETSc and MOAB data structures
1101 const char param_file[] = "param_file.petsc";
1102 MoFEM::Core::Initialize(&argc, &argv, param_file, help);
1103
1104 // Add logging channel for example
1105 auto core_log = logging::core::get();
1106 core_log->add_sink(
1108 LogManager::setLog("Seepage");
1109 MOFEM_LOG_TAG("Seepage", "seepage");
1110 core_log->add_sink(
1112 LogManager::setLog("SeepageSync");
1113 MOFEM_LOG_TAG("SeepageSync", "seepage");
1114
1115 try {
1116
1117 //! [Register MoFEM discrete manager in PETSc]
1118 DMType dm_name = "DMMOFEM";
1119 CHKERR DMRegister_MoFEM(dm_name);
1120 //! [Register MoFEM discrete manager in PETSc
1121
1122 //! [Create MoAB]
1123 moab::Core mb_instance; ///< mesh database
1124 moab::Interface &moab = mb_instance; ///< mesh database interface
1125 //! [Create MoAB]
1126
1127 //! [Create MoFEM]
1128 MoFEM::Core core(moab); ///< finite element database
1129 MoFEM::Interface &m_field = core; ///< finite element database insterface
1130 //! [Create MoFEM]
1131
1132 //! [Load mesh]
1133 Simple *simple = m_field.getInterface<Simple>();
1135 CHKERR simple->loadFile();
1136 //! [Load mesh]
1137
1138 //! [Seepage]
1139 Seepage ex(m_field);
1140 CHKERR ex.runProblem();
1141 //! [Seepage]
1142 }
1144
1146}
static auto filter_true_skin(MoFEM::Interface &m_field, Range &&skin)
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
int main()
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
[Define dimension]
Kronecker Delta class symmetric.
@ QUIET
@ ROW
#define CATCH_ERRORS
Catch errors.
@ MF_EXIST
FieldApproximationBase
approximation base
Definition definitions.h:58
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
Definition definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ L2
field with C-1 continuity
Definition definitions.h:88
@ H1
continuous field
Definition definitions.h:85
@ NOSPACE
Definition definitions.h:83
@ HCURL
field with continuous tangents
Definition definitions.h:86
@ HDIV
field with continuous normal traction
Definition definitions.h:87
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ TEMPERATURESET
@ HEATFLUXSET
@ NODESET
@ SIDESET
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
constexpr int order
double bulk_modulus_K
double shear_modulus_G
auto integration_rule
constexpr auto t_kd
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMixScalarTimesDiv< SPACE_DIM, COORD_TYPE > OpMixScalarTimesDiv
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
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:576
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1102
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition DMMoFEM.hpp:1059
SmartPetscObj< TS > createTSIM(SmartPetscObj< DM > dm=nullptr)
Create TS (time) implicit solver.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
SeverityLevel
Severity levels.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
FTensor::Index< 'i', SPACE_DIM > i
PetscBool doEvalField
const double c
speed of light (cm/ns)
double D
const double n
refractive index of diffusive medium
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
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:1046
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition DMMoFEM.hpp:1144
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *ctx)
Sens monitor printing residual field by field.
Definition SnesCtx.cpp:592
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
static FTensor::Ddg< FTensor::PackPtr< T *, 1 >, Tensor_Dim01, Tensor_Dim23 > getFTensor4DdgFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 4 on first two and last indices from form data matrix.
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition DMMoFEM.hpp:1130
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)
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
double young_modulus
Young modulus.
Definition plastic.cpp:125
double poisson_ratio
Poisson ratio.
Definition plastic.cpp:126
constexpr auto size_symm
Definition plastic.cpp:42
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition radiation.cpp:29
auto save_range
Definition seepage.cpp:169
double default_conductivity
Definition seepage.cpp:163
OpBaseDotH OpBaseDivFlux
Integrate Rhs base of temperature times divergent of flux (T)
Definition seepage.cpp:129
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition seepage.cpp:50
static char help[]
[Solve]
Definition seepage.cpp:1096
double fluid_density
Definition seepage.cpp:164
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpBaseDotH
Integrate Rhs base of temperature time heat capacity times heat rate (T)
Definition seepage.cpp:122
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:48
#define EXECUTABLE_DIMENSION
Definition seepage.cpp:24
constexpr int SPACE_DIM
Definition seepage.cpp:33
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixScalarTimesDiv< SPACE_DIM > OpBaseDivU
[Essential boundary conditions]
Definition seepage.cpp:77
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 3, 3, 1 > OpHdivFlux
Integrating Rhs flux base (1/k) flux (FLUX)
Definition seepage.cpp:107
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:85
double scale
[Essential boundary conditions (Least square approach)]
Definition seepage.cpp:159
double default_young_modulus
Definition seepage.cpp:161
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:93
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition seepage.cpp:62
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 1, 2 > OpHDivH
Integrate Rhs div flux base times temperature (T)
Definition seepage.cpp:114
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition seepage.cpp:64
double default_poisson_ratio
Definition seepage.cpp:162
FTensor::Index< 'm', 3 > m
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition BcManager.hpp:29
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Core (interface) class.
Definition Core.hpp:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:118
Deprecated interface functions.
Definition of the displacement bc data structure.
Definition BCData.hpp:72
Data on single entity (This is passed as argument to DataOperator::doWork)
Essential boundary conditions.
Class (Function) to enforce essential constrains on the left hand side diagonal.
Definition Essential.hpp:33
Class (Function) to enforce essential constrains on the right hand side diagonal.
Definition Essential.hpp:41
Class (Function) to calculate residual side diagonal.
Definition Essential.hpp:49
Class (Function) to enforce essential constrains.
Definition Essential.hpp:25
structure for User Loop Methods on finite elements
Field evaluator interface.
Definition of the heat flux bc data structure.
Definition BCData.hpp:423
Section manager is used to create indexes and sections.
Definition ISManager.hpp:23
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.
Interface for managing meshsets containing materials and boundary conditions.
Assembly methods.
Definition Natural.hpp:65
Get vector field for H-div approximation.
Calculate divergence of vector field.
Get value at integration points for scalar field.
Calculates the trace of an input matrix.
Get field gradients time derivative at integration pts for scalar field rank 0, i....
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Get values at integration pts for tensor field rank 1, i.e. vector field.
Enforce essential constrains on lhs.
Definition Essential.hpp:81
Enforce essential constrains on rhs.
Definition Essential.hpp:65
Post post-proc data at points from hash maps.
Set indices on entities on finite element.
PipelineManager interface.
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition Simple.hpp:27
MoFEMErrorCode addDomainField(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition Simple.cpp:261
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Seepage(MoFEM::Interface &m_field)
Definition seepage.cpp:180
MoFEMErrorCode setupProblem()
[Run problem]
Definition seepage.cpp:448
std::array< double, SPACE_DIM > fieldEvalCoords
Definition seepage.cpp:194
MoFEMErrorCode createCommonData()
[Set up problem]
Definition seepage.cpp:475
MoFEMErrorCode bC()
[Create common data]
Definition seepage.cpp:512
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition seepage.cpp:198
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition seepage.cpp:196
MoFEMErrorCode runProblem()
[Run problem]
Definition seepage.cpp:436
MoFEMErrorCode OPs()
[Boundary condition]
Definition seepage.cpp:613
PetscBool doEvalField
Definition seepage.cpp:193
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition seepage.cpp:197
MoFEM::Interface & mField
Definition seepage.cpp:185
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:220
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
Definition seepage.cpp:792