v0.14.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>
31
32using namespace MoFEM;
33
34constexpr int SPACE_DIM =
35 EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
36
44
47
48//! [Only used with Hooke equation (linear material model)]
50 GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
53//! [Only used with Hooke equation (linear material model)]
54
55//! [Only used for dynamics]
60//! [Only used for dynamics]
61
62//! [Only used with Hencky/nonlinear material]
64 GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
67//! [Only used with Hencky/nonlinear material]
68
69//! [Essential boundary conditions]
76//! [Essential boundary conditions]
78
81
82// Thermal operators
83/**
84 * @brief Integrate Lhs base of flux (1/k) base of flux (FLUX x FLUX)
85 *
86 */
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 */
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 */
105
106/**
107 * @brief Integrating Rhs flux base (1/k) flux (FLUX)
108 */
110 GAUSS>::OpBaseTimesVector<3, 3, 1>;
111
112/**
113 * @brief Integrate Rhs div flux base times temperature (T)
114 *
115 */
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 */
125 GAUSS>::OpBaseTimesScalarField<1>;
126
127/**
128 * @brief Integrate Rhs base of temerature times divergenc of flux (T)
129 *
130 */
132
133//! [Body and heat source]
142//! [Body and heat source]
143
144//! [Natural boundary conditions]
150//! [Natural boundary conditions]
151
152//! [Essential boundary conditions (Least square approach)]
155 GAUSS>::OpEssentialRhs<HeatFluxCubitBcData, 3, SPACE_DIM>;
158 GAUSS>::OpEssentialLhs<HeatFluxCubitBcData, 3, SPACE_DIM>;
159//! [Essential boundary conditions (Least square approach)]
160
161double scale = 1.;
162
164double default_poisson_ratio = 0.25; // zero for verification
166double fluid_density = 10;
167
168// #include <OpPostProcElastic.hpp>
169#include <SeepageOps.hpp>
170
171auto save_range = [](moab::Interface &moab, const std::string name,
172 const Range r) {
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);
178};
179
180struct Seepage {
181
182 Seepage(MoFEM::Interface &m_field) : mField(m_field) {}
183
185
186private:
188
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
203 : public boost::enable_shared_from_this<BlockedParameters> {
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
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
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) {
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)
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,
243
244 for (auto &b : blockData) {
245
246 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
247 CHKERR getMatDPtr(matDPtr, b.bulkModulusK, b.shearModulusG);
249 }
250 }
251
252 CHKERR getMatDPtr(matDPtr, bulkModulusKDefault, shearModulusGDefault);
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
270 extractElasticBlockData(MoFEM::Interface &m_field,
271 std::vector<const CubitMeshSets *> meshset_vec_ptr,
272 Sev sev) {
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");
304 }
305
306 MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
307 double bulk_modulus_K, double shear_modulus_G) {
309 //! [Calculate elasticity tensor]
310 auto set_material_stiffness = [&]() {
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();
331 }
332 };
333
334 double default_bulk_modulus_K =
336 double default_shear_modulus_G =
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)
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,
365
366 for (auto &b : blockData) {
367
368 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
369 *conductivityPtr = b.conductivity;
371 }
372 }
373
374 *conductivityPtr = default_conductivity;
375
377 }
378
379 private:
380 struct BlockData {
381 double conductivity;
382 Range blockEnts;
383 };
384
385 std::vector<BlockData> blockData;
386
388 extractThermallockData(MoFEM::Interface &m_field,
389 std::vector<const CubitMeshSets *> meshset_vec_ptr,
390 Sev sev) {
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");
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
427
428 (boost::format("%s(.*)") % block_thermal_name).str()
429
430 ))
431
432 ));
433
435}
436
437//! [Run problem]
442 CHKERR bC();
443 CHKERR OPs();
444 CHKERR tsSolve();
446}
447//! [Run problem]
448
449//! [Set up problem]
453 // Add field
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
473}
474//! [Set up problem]
475
476//! [Create common data]
479
480 auto get_command_line_parameters = [&]() {
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
505 };
506
507 CHKERR get_command_line_parameters();
508
510}
511//! [Create common data]
512
513//! [Boundary condition]
516
517 MOFEM_LOG("SYNC", Sev::noisy) << "bC";
519
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) {
540 for (auto m :
541
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 }
551 };
552
553 auto remove_named_blocks = [&](auto n) {
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 }
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;
597
598#ifndef NDEBUG
599
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
611}
612//! [Boundary condition]
613
614//! [Push operators to pipeline]
617 auto pipeline_mng = mField.getInterface<PipelineManager>();
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) {
651
652 CHKERR addMatBlockOps(pip, "MAT_ELASTIC", "MAT_FLUID", block_params,
653 Sev::inform);
655
657 "U", u_grad_ptr));
658 pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(u_grad_ptr, strain_ptr));
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));
668 "FLUX", div_flux_ptr));
669
671 };
672
673 auto add_domain_ops_lhs_mechanical = [&](auto &pip) {
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));
681 };
682
683 auto add_domain_ops_rhs_mechanical = [&](auto &pip) {
685
687 pip, mField, "U", {time_scale}, "BODY_FORCE", Sev::inform);
688
689 // Calculate internal forece
691 "U", strain_ptr, stress_ptr, mDPtr));
692 pip.push_back(new OpInternalForceCauchy("U", stress_ptr));
693 pip.push_back(
695
697 };
698
699 auto add_domain_ops_lhs_seepage = [&](auto &pip, auto &fe) {
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
718 };
719
720 auto add_domain_ops_rhs_seepage = [&](auto &pip) {
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
733 };
734
735 auto add_boundary_rhs_ops = [&](auto &pip) {
737
739
740 pip.push_back(new OpSetBc("FLUX", true, boundary_marker));
741
743 pipeline_mng->getOpBoundaryRhsPipeline(), mField, "U", {time_scale},
744 "FORCE", Sev::inform);
745
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));
757 mField, pip, simple->getProblemName(), "FLUX", mat_flux_ptr,
758 {time_scale});
759
761 };
762
763 auto add_boundary_lhs_ops = [&](auto &pip) {
765
767
770 mField, pip, simple->getProblemName(), "FLUX");
771
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
790}
791//! [Push operators to pipeline]
792
793//! [Solve]
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) {
806 SNES snes;
807 CHKERR TSGetSNES(solver, &snes);
808 CHKERR SNESMonitorSet(snes,
809 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
811 (void *)(snes_ctx_ptr.get()), nullptr);
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);
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(
840 post_proc_fe->getOpPtrVector().push_back(
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(
847 "U", mat_strain_ptr, mat_stress_ptr, mDPtr));
848
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);
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
889 "U", u_grad_ptr));
890 pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(u_grad_ptr, strain_ptr));
891
892 // Calculate internal forece
894 "U", strain_ptr, stress_ptr, mDPtr));
895 pip.push_back(new OpInternalForceCauchy("U", stress_ptr));
896
897 fe_ptr->postProcessHook =
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) {
910 monitor_ptr->preProcessHook = [&]() {
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);
964 field_eval_ptr->getOpPtrVector(), {H1, HDIV});
965 field_eval_ptr->getOpPtrVector().push_back(
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(
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
1000 };
1001 auto null = boost::shared_ptr<FEMethod>();
1002 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(), null,
1003 monitor_ptr, null);
1005 };
1006
1007 auto set_fieldsplit_preconditioner = [&](auto solver) {
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
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]() {
1074 mField, post_proc_rhs_ptr, nullptr, Sev::verbose)();
1076 mField, post_proc_rhs_ptr, 1.)();
1078 };
1079 auto get_post_proc_hook_lhs = [this, post_proc_lhs_ptr]() {
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
1106}
1107//! [Solve]
1108
1109static char help[] = "...\n\n";
1110
1111int 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(
1121 LogManager::setLog("Seepage");
1122 MOFEM_LOG_TAG("Seepage", "seepage");
1123 core_log->add_sink(
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 }
1157
1159}
static Index< 'o', 3 > o
std::string param_file
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
Definition: LogManager.hpp:352
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
int main()
Definition: adol-c_atom.cpp:46
constexpr int SPACE_DIM
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
@ QUIET
Definition: definitions.h:208
@ ROW
Definition: definitions.h:123
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ MF_EXIST
Definition: definitions.h:100
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ 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
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ TEMPERATURESET
Definition: definitions.h:155
@ HEATFLUXSET
Definition: definitions.h:156
@ NODESET
Definition: definitions.h:146
@ SIDESET
Definition: definitions.h:147
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
constexpr int order
double bulk_modulus_K
double shear_modulus_G
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
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:47
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:572
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1003
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
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
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
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:1042
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1045
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
Definition: SnesCtx.cpp:232
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)
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition: Templates.hpp:1857
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition: DMMoFEM.hpp:1031
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
constexpr AssemblyType A
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
double young_modulus
Young modulus.
Definition: plastic.cpp:172
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:171
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: seepage.cpp:66
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: seepage.cpp:64
double default_conductivity
Definition: seepage.cpp:165
OpBaseDotH OpBaseDivFlux
Integrate Rhs base of temerature times divergenc of flux (T)
Definition: seepage.cpp:131
static char help[]
[Solve]
Definition: seepage.cpp:1109
double fluid_density
Definition: seepage.cpp:166
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: seepage.cpp:52
#define EXECUTABLE_DIMENSION
Definition: seepage.cpp:24
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixScalarTimesDiv< SPACE_DIM > OpBaseDivU
Definition: seepage.cpp:80
constexpr int SPACE_DIM
Definition: seepage.cpp:34
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 3, 3, 1 > OpHdivFlux
Integrating Rhs flux base (1/k) flux (FLUX)
Definition: seepage.cpp:110
double scale
[Essential boundary conditions (Least square approach)]
Definition: seepage.cpp:161
double default_young_modulus
Definition: seepage.cpp:163
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
Definition: seepage.cpp:73
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
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
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 1, 2 > OpHDivH
Integrate Rhs div flux base times temperature (T)
Definition: seepage.cpp:117
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
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
Definition: seepage.cpp:75
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: seepage.cpp:71
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
double default_poisson_ratio
Definition: seepage.cpp:164
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
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:112
Deprecated interface functions.
Definition of the displacement bc data structure.
Definition: BCData.hpp:76
Data on single entity (This is passed as argument to DataOperator::doWork)
Essential boundary conditions.
Definition: Essential.hpp:111
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 enforce essential constrains.
Definition: Essential.hpp:25
Class (Function) to calculate residual side diagonal.
Definition: Essential.hpp:49
structure for User Loop Methods on finite elements
Field evaluator interface.
@ OPSPACE
operator do Work is execute on space data
Definition of the heat flux bc data structure.
Definition: BCData.hpp:427
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.
Definition: LogManager.cpp:298
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.
Definition: LogManager.cpp:348
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 filed rank 0, i....
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed 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.
Scale base functions by inverses of measure of element.
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
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Seepage(MoFEM::Interface &m_field)
Definition: seepage.cpp:182
MoFEMErrorCode setupProblem()
[Run problem]
Definition: seepage.cpp:450
std::array< double, SPACE_DIM > fieldEvalCoords
Definition: seepage.cpp:196
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: seepage.cpp:477
MoFEMErrorCode bC()
[Create common data]
Definition: seepage.cpp:514
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: seepage.cpp:200
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: seepage.cpp:198
MoFEMErrorCode runProblem()
[Run problem]
Definition: seepage.cpp:438
MoFEMErrorCode OPs()
[Boundary condition]
Definition: seepage.cpp:615
PetscBool doEvalField
Definition: seepage.cpp:195
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: seepage.cpp:199
MoFEM::Interface & mField
Definition: seepage.cpp:187
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
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
Definition: seepage.cpp:794
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 3, SPACE_DIM, 1 > OpHdivFlux
Integrating Rhs flux base (1/k) flux (FLUX)
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)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, SPACE_DIM > OpHdivHdiv
[Linear elastic problem]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Linear elastic problem]
OpBaseDotT OpBaseDivFlux
Integrate Rhs base of temperature times divergence of flux (T)