13#include <petsc/private/petscimpl.h>
17static char help[] =
"...\n\n";
19#undef PYTHON_INIT_SURFACE
21#ifdef PYTHON_INIT_SURFACE
22#include <boost/python.hpp>
23#include <boost/python/def.hpp>
24namespace bp = boost::python;
27 SurfacePython() =
default;
28 virtual ~SurfacePython() =
default;
46 auto main_module = bp::import(
"__main__");
47 mainNamespace = main_module.attr(
"__dict__");
48 bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
50 surfaceFun = mainNamespace[
"surface"];
52 }
catch (bp::error_already_set
const &) {
75 double x,
double y,
double z,
double eta,
double &s
82 s = bp::extract<double>(surfaceFun(x, y, z,
eta));
84 }
catch (bp::error_already_set
const &) {
93 bp::object mainNamespace;
94 bp::object surfaceFun;
97static boost::weak_ptr<SurfacePython> surfacePythonWeakPtr;
109 IntegrationType::GAUSS;
149template <CoordinateTypes COORD_TYPE>
202double tol = std::numeric_limits<float>::epsilon();
237 constexpr int sub_stepping = 16;
238 return std::min(1.,
static_cast<double>(ts_step) / sub_stepping);
248auto my_max = [](
const double x) {
return (x - 1 + std::abs(x + 1)) / 2; };
255auto my_min = [](
const double x) {
return (x + 1 - std::abs(x - 1)) / 2; };
272 if (
h >= -1 &&
h < 1)
289 return diff *
h + ave;
312auto get_f = [](
const double h) {
return 4 *
W *
h * (
h *
h - 1); };
343 return md * (1 -
h *
h);
361 const double h2 =
h *
h;
362 const double h3 = h2 *
h;
364 return md * (2 * h3 - 3 * h2 + 1);
366 return md * (-2 * h3 - 3 * h2 + 1);
376 return md * (6 *
h * (
h - 1));
378 return md * (-6 *
h * (
h + 1));
416 constexpr double R = 0.0125;
417 constexpr double A =
R * 0.2;
418 const double theta = atan2(r, y);
419 const double w =
R +
A * cos(
n * theta);
420 const double d = std::sqrt(r * r + y * y);
421 return tanh((w - d) / (
eta * std::sqrt(2)));
434 constexpr double y0 = 0.5;
435 constexpr double R = 0.4;
437 const double d = std::sqrt(r * r + y * y);
438 return tanh((
R - d) / (
eta * std::sqrt(2)));
451 constexpr double water_height = 0.;
452 return tanh((water_height - y) / (
eta * std::sqrt(2)));
466 return -tanh((-0.039 - x) / (
eta * std::sqrt(2)));
476auto init_h = [](
double r,
double y,
double theta) {
477#ifdef PYTHON_INIT_SURFACE
479 if (
auto ptr = surfacePythonWeakPtr.lock()) {
481 "error eval python");
521 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel();
531 MPI_Allreduce(&l_size, &g_size, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
544auto save_range = [](moab::Interface &moab,
const std::string name,
549 CHKERR moab.add_entities(*out_meshset, r);
550 CHKERR moab.write_file(name.c_str(),
"MOAB",
"PARALLEL=WRITE_PART",
551 out_meshset->get_ptr(), 1);
568 std::vector<EntityHandle> ents_vec;
574 auto dofs = prb_ptr->numeredRowDofsPtr;
577 ents_vec.reserve(std::distance(lo_it, hi_it));
579 for (; lo_it != hi_it; ++lo_it) {
580 ents_vec.push_back((*lo_it)->getEnt());
583 std::sort(ents_vec.begin(), ents_vec.end());
584 auto it = std::unique(ents_vec.begin(), ents_vec.end());
586 r.insert_list(ents_vec.begin(), it);
600 std::vector<EntityHandle> ents_vec;
605 auto dofs = prb_ptr->numeredRowDofsPtr;
606 ents_vec.reserve(dofs->size());
608 for (
auto d : *dofs) {
609 ents_vec.push_back(d->getEnt());
612 std::sort(ents_vec.begin(), ents_vec.end());
613 auto it = std::unique(ents_vec.begin(), ents_vec.end());
615 r.insert_list(ents_vec.begin(), it);
744 PetscReal
a, Mat
A, Mat
B,
788 boost::shared_ptr<SnesCtx>
790 boost::shared_ptr<TsCtx>
949 boost::shared_ptr<FEMethod> fe_top, std::string
field_name,
950 ForcesAndSourcesCore::UserDataOperator::OpType op,
952 boost::function<boost::shared_ptr<ForcesAndSourcesCore>()> get_elem,
973 MOFEM_LOG(
"FS", Sev::inform) <<
"Read mesh for problem";
996 const char *coord_type_names[] = {
"cartesian",
"polar",
"cylindrical",
1001 MOFEM_LOG(
"FS", Sev::inform) <<
"Approximation order = " <<
order;
1003 <<
"Number of refinement levels nb_levels = " <<
nb_levels;
1010 "-rho_m", &
rho_m, PETSC_NULLPTR);
1012 &
mu_m, PETSC_NULLPTR);
1014 &
rho_p, PETSC_NULLPTR);
1016 &
mu_p, PETSC_NULLPTR);
1020 "Height of the well in energy functional",
"-W",
1034 MOFEM_LOG(
"FS", Sev::inform) <<
"Acceleration a0 = " <<
a0;
1035 MOFEM_LOG(
"FS", Sev::inform) <<
"\"Minus\" phase density rho_m = " <<
rho_m;
1036 MOFEM_LOG(
"FS", Sev::inform) <<
"\"Minus\" phase viscosity mu_m = " <<
mu_m;
1037 MOFEM_LOG(
"FS", Sev::inform) <<
"\"Plus\" phase density rho_p = " <<
rho_p;
1038 MOFEM_LOG(
"FS", Sev::inform) <<
"\"Plus\" phase viscosity mu_p = " <<
mu_p;
1039 MOFEM_LOG(
"FS", Sev::inform) <<
"Surface tension lambda = " <<
lambda;
1041 <<
"Height of the well in energy functional W = " <<
W;
1042 MOFEM_LOG(
"FS", Sev::inform) <<
"Characteristic mesh size h = " <<
h;
1043 MOFEM_LOG(
"FS", Sev::inform) <<
"Mobility md = " <<
md;
1047 MOFEM_LOG(
"FS", Sev::inform) <<
"Average viscosity mu_ave = " <<
mu_ave;
1048 MOFEM_LOG(
"FS", Sev::inform) <<
"Difference viscosity mu_diff = " <<
mu_diff;
1078 auto set_problem_bit = [&]() {
1089 CHKERR set_problem_bit();
1101#ifdef PYTHON_INIT_SURFACE
1102 auto get_py_surface_init = []() {
1103 auto py_surf_init = boost::make_shared<SurfacePython>();
1104 CHKERR py_surf_init->surfaceInit(
"surface.py");
1105 surfacePythonWeakPtr = py_surf_init;
1106 return py_surf_init;
1108 auto py_surf_init = get_py_surface_init();
1115 auto dm =
simple->getDM();
1119 auto reset_bits = [&]() {
1123 start_mask[s] =
true;
1125 CHKERR bit_mng->lambdaBitRefLevel(
1134 auto add_parent_field = [&](
auto fe,
auto op,
auto field) {
1139 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1147 auto h_ptr = boost::make_shared<VectorDouble>();
1148 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
1149 auto g_ptr = boost::make_shared<VectorDouble>();
1150 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
1152 auto set_generic = [&](
auto fe) {
1154 auto &pip = fe->getOpPtrVector();
1162 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1165 fe_parent->getOpPtrVector(), {H1});
1171 CHKERR add_parent_field(fe, UDO::OPROW,
"H");
1176 CHKERR add_parent_field(fe, UDO::OPROW,
"G");
1184 auto post_proc = [&](
auto exe_test) {
1186 auto post_proc_fe = boost::make_shared<PostProcEleDomain>(
mField);
1187 post_proc_fe->exeTestHook = exe_test;
1189 CHKERR set_generic(post_proc_fe);
1193 post_proc_fe->getOpPtrVector().push_back(
1195 new OpPPMap(post_proc_fe->getPostProcMesh(),
1196 post_proc_fe->getMapGaussPts(),
1198 {{
"H", h_ptr}, {
"G", g_ptr}},
1200 {{
"GRAD_H", grad_h_ptr}, {
"GRAD_G", grad_g_ptr}},
1211 CHKERR post_proc_fe->writeFile(
"out_init.h5m");
1216 auto solve_init = [&](
auto exe_test) {
1219 pip_mng->getOpDomainRhsPipeline().clear();
1220 pip_mng->getOpDomainLhsPipeline().clear();
1222 auto set_domain_rhs = [&](
auto fe) {
1225 auto &pip = fe->getOpPtrVector();
1227 CHKERR add_parent_field(fe, UDO::OPROW,
"H");
1228 pip.push_back(
new OpRhsH<true>(
"H",
nullptr,
nullptr, h_ptr, grad_h_ptr,
1230 CHKERR add_parent_field(fe, UDO::OPROW,
"G");
1231 pip.push_back(
new OpRhsG<true>(
"G", h_ptr, grad_h_ptr, g_ptr));
1235 auto set_domain_lhs = [&](
auto fe) {
1239 auto &pip = fe->getOpPtrVector();
1241 CHKERR add_parent_field(fe, UDO::OPROW,
"H");
1242 CHKERR add_parent_field(fe, UDO::OPCOL,
"H");
1245 CHKERR add_parent_field(fe, UDO::OPCOL,
"G");
1248 CHKERR add_parent_field(fe, UDO::OPROW,
"G");
1251 CHKERR add_parent_field(fe, UDO::OPCOL,
"H");
1256 auto create_subdm = [&]() {
1257 auto level_ents_ptr = boost::make_shared<Range>();
1263 CHKERR DMSetType(subdm,
"DMMOFEM");
1269 for (
auto f : {
"H",
"G"}) {
1275 if constexpr (
debug) {
1279 dm_ents.subset_by_type(MBVERTEX));
1280 dm_ents = subtract(dm_ents, dm_ents.subset_by_type(MBVERTEX));
1288 auto subdm = create_subdm();
1290 pip_mng->getDomainRhsFE().reset();
1291 pip_mng->getDomainLhsFE().reset();
1294 pip_mng->getDomainLhsFE()->exeTestHook = exe_test;
1295 pip_mng->getDomainRhsFE()->exeTestHook = exe_test;
1297 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
1298 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
1301 auto snes = pip_mng->createSNES(subdm);
1304 auto set_section_monitor = [&](
auto solver) {
1306 CHKERR SNESMonitorSet(solver,
1309 (
void *)(snes_ctx_ptr.get()),
nullptr);
1314 auto print_section_field = [&]() {
1319 PetscInt num_fields;
1320 CHKERR PetscSectionGetNumFields(section, &num_fields);
1321 for (
int f = 0;
f < num_fields; ++
f) {
1325 <<
"Field " <<
f <<
" " << std::string(
field_name);
1331 CHKERR set_section_monitor(snes);
1332 CHKERR print_section_field();
1334 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
1338 CHKERR SNESSetFromOptions(snes);
1339 CHKERR SNESSolve(snes, PETSC_NULLPTR,
D);
1341 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1342 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1353 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
1364 pip_mng->getOpDomainRhsPipeline().clear();
1365 pip_mng->getOpDomainLhsPipeline().clear();
1368 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM_X",
1370 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM_X",
1372 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM_Y",
1374 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM_Y",
1376 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX",
"U",
1378 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX",
"L",
1380 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"ZERO",
1402 auto fe_domain_lhs = pip_mng->getDomainLhsFE();
1403 auto fe_domain_rhs = pip_mng->getDomainRhsFE();
1404 auto fe_bdy_lhs = pip_mng->getBoundaryLhsFE();
1405 auto fe_bdy_rhs = pip_mng->getBoundaryRhsFE();
1407 pip_mng->getDomainLhsFE().reset();
1408 pip_mng->getDomainRhsFE().reset();
1409 pip_mng->getBoundaryLhsFE().reset();
1410 pip_mng->getBoundaryRhsFE().reset();
1415 auto add_parent_field_domain = [&](
auto fe,
auto op,
auto field,
auto bit) {
1420 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1429 auto add_parent_field_bdy = [&](
auto fe,
auto op,
auto field,
auto bit) {
1434 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1444 auto create_run_parent_op = [&](
auto parent_fe_ptr,
auto this_fe_ptr,
1446 auto parent_mask = fe_bit;
1454 auto get_parents_vel_fe_ptr = [&](
auto this_fe_ptr,
auto fe_bit) {
1455 std::vector<boost::shared_ptr<DomainParentEle>> parents_elems_ptr_vec;
1457 parents_elems_ptr_vec.emplace_back(
1458 boost::make_shared<DomainParentEle>(
mField));
1460 parents_elems_ptr_vec[
l - 1]->getOpPtrVector().push_back(
1461 create_run_parent_op(parents_elems_ptr_vec[
l], this_fe_ptr, fe_bit));
1463 return parents_elems_ptr_vec[0];
1467 auto solve_projection = [&](
auto exe_test) {
1470 auto set_domain_rhs = [&](
auto fe) {
1472 auto &pip = fe->getOpPtrVector();
1474 auto u_ptr = boost::make_shared<MatrixDouble>();
1475 auto p_ptr = boost::make_shared<VectorDouble>();
1476 auto h_ptr = boost::make_shared<VectorDouble>();
1477 auto g_ptr = boost::make_shared<VectorDouble>();
1479 auto eval_fe_ptr = boost::make_shared<DomainParentEle>(
mField);
1481 auto &pip = eval_fe_ptr->getOpPtrVector();
1487 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1495 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW,
"U",
1498 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW,
"P",
1501 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW,
"H",
1504 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW,
"G",
1508 auto parent_eval_fe_ptr =
1510 pip.push_back(create_run_parent_op(parent_eval_fe_ptr, eval_fe_ptr,
1513 auto assemble_fe_ptr = boost::make_shared<DomainParentEle>(
mField);
1515 auto &pip = assemble_fe_ptr->getOpPtrVector();
1521 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1527 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW,
"U",
1530 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW,
"P",
1533 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW,
"H",
1536 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW,
"G",
1540 auto parent_assemble_fe_ptr =
1542 pip.push_back(create_run_parent_op(
1548 auto set_domain_lhs = [&](
auto fe) {
1551 auto &pip = fe->getOpPtrVector();
1559 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1568 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"U",
1570 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"U",
1573 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"P",
1575 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"P",
1578 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"H",
1580 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"H",
1583 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"G",
1585 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"G",
1592 auto set_bdy_rhs = [&](
auto fe) {
1594 auto &pip = fe->getOpPtrVector();
1596 auto l_ptr = boost::make_shared<VectorDouble>();
1598 auto eval_fe_ptr = boost::make_shared<BoundaryParentEle>(
mField);
1600 auto &pip = eval_fe_ptr->getOpPtrVector();
1605 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1613 CHKERR add_parent_field_bdy(eval_fe_ptr, UDO::OPROW,
"L",
1617 auto parent_eval_fe_ptr =
1619 pip.push_back(create_run_parent_op(parent_eval_fe_ptr, eval_fe_ptr,
1622 auto assemble_fe_ptr = boost::make_shared<BoundaryParentEle>(
mField);
1624 auto &pip = assemble_fe_ptr->getOpPtrVector();
1629 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1637 OpLSize(boost::shared_ptr<VectorDouble> l_ptr)
1641 if (lPtr->size() != getGaussPts().size2()) {
1642 lPtr->resize(getGaussPts().size2());
1649 boost::shared_ptr<VectorDouble> lPtr;
1652 pip.push_back(
new OpLSize(l_ptr));
1654 CHKERR add_parent_field_bdy(assemble_fe_ptr, UDO::OPROW,
"L",
1658 auto parent_assemble_fe_ptr =
1660 pip.push_back(create_run_parent_op(
1666 auto set_bdy_lhs = [&](
auto fe) {
1669 auto &pip = fe->getOpPtrVector();
1675 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1684 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"L",
1686 CHKERR add_parent_field_bdy(fe, UDO::OPCOL,
"L",
1694 auto level_ents_ptr = boost::make_shared<Range>();
1698 auto get_prj_ents = [&]() {
1703 auto common_ents = intersect(prj_mesh, *level_ents_ptr);
1704 prj_mesh = subtract(unite(*level_ents_ptr, prj_mesh), common_ents)
1710 auto prj_ents = get_prj_ents();
1714 auto rebuild = [&]() {
1718 std::vector<std::string> fields{
"U",
"P",
"H",
"G",
"L"};
1719 std::map<std::string, boost::shared_ptr<Range>> range_maps{
1721 {
"U", level_ents_ptr},
1722 {
"P", level_ents_ptr},
1723 {
"H", level_ents_ptr},
1724 {
"G", level_ents_ptr},
1725 {
"L", level_ents_ptr}
1730 simple->getProblemName(), PETSC_TRUE,
1731 &range_maps, &range_maps);
1734 CHKERR prb_mng->partitionFiniteElements(
"SUB_SOLVER",
true, 0,
1737 CHKERR prb_mng->partitionGhostDofsOnDistributedMesh(
"SUB_SOLVER");
1742 MOFEM_LOG(
"FS", Sev::verbose) <<
"Create projection problem";
1746 auto dm =
simple->getDM();
1749 CHKERR DMSetType(subdm,
"DMMOFEM");
1754 MOFEM_LOG(
"FS", Sev::inform) <<
"Nothing to project";
1759 auto create_dummy_dm = [&]() {
1762 simple->getProblemName().c_str(),
1768 auto subdm = create_subdm();
1771 pip_mng->getDomainRhsFE().reset();
1772 pip_mng->getDomainLhsFE().reset();
1773 pip_mng->getBoundaryRhsFE().reset();
1774 pip_mng->getBoundaryLhsFE().reset();
1779 pip_mng->getDomainLhsFE()->exeTestHook = exe_test;
1780 pip_mng->getDomainRhsFE()->exeTestHook = [](
FEMethod *fe_ptr) {
1783 pip_mng->getBoundaryLhsFE()->exeTestHook = exe_test;
1784 pip_mng->getBoundaryRhsFE()->exeTestHook = [](
FEMethod *fe_ptr) {
1788 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
1789 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
1790 CHKERR set_bdy_rhs(pip_mng->getCastBoundaryRhsFE());
1791 CHKERR set_bdy_lhs(pip_mng->getCastBoundaryLhsFE());
1796 auto ksp = pip_mng->createKSP(subdm);
1797 CHKERR KSPSetFromOptions(ksp);
1801 auto get_norm = [&](
auto x) {
1803 CHKERR VecNorm(x, NORM_2, &nrm);
1810 auto zero_dofs = [](boost::shared_ptr<FieldEntity> ent_ptr) {
1812 for (
auto &
v : ent_ptr->getEntFieldData()) {
1818 auto solve = [&](
auto S) {
1820 CHKERR VecZeroEntries(S);
1822 CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
1823 CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
1825 CHKERR VecGhostUpdateBegin(S, INSERT_VALUES, SCATTER_FORWARD);
1826 CHKERR VecGhostUpdateEnd(S, INSERT_VALUES, SCATTER_FORWARD);
1833 MOFEM_LOG(
"FS", Sev::inform) <<
"Solve projection";
1838 auto dummy_dm = create_dummy_dm();
1843 auto apply_restrict = [&]() {
1844 auto get_is = [](
auto v) {
1846 auto create = [&]() {
1850 CHKERR VecGetOwnershipRange(
v, &ystart, NULL);
1851 CHKERR ISCreateStride(PETSC_COMM_SELF,
n, ystart, 1, &iy);
1858 auto iy = get_is(glob_x);
1862 DMSubDomainRestrict(
simple->getDM(), s, PETSC_NULLPTR, dummy_dm),
1868 DMGetNamedGlobalVector(dummy_dm,
"TSTheta_Xdot", &Xdot),
1871 auto forward_ghost = [](
auto g) {
1873 CHKERR VecGhostUpdateBegin(
g, INSERT_VALUES, SCATTER_FORWARD);
1874 CHKERR VecGhostUpdateEnd(
g, INSERT_VALUES, SCATTER_FORWARD);
1881 if constexpr (
debug) {
1883 <<
"Reverse restrict: X0 " << get_norm(X0) <<
" Xdot "
1887 return std::vector<Vec>{X0, Xdot};
1890 auto ts_solver_vecs = apply_restrict();
1892 if (ts_solver_vecs.size()) {
1894 for (
auto v : ts_solver_vecs) {
1895 MOFEM_LOG(
"FS", Sev::inform) <<
"Solve projection vector";
1901 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
1902 MOFEM_LOG(
"WORLD", Sev::verbose) <<
"Zero field " << f;
1903 CHKERR field_blas->fieldLambdaOnEntities(zero_dofs, f);
1911 MOFEM_LOG(
"FS", Sev::inform) <<
"Norm V " << get_norm(
v);
1914 CHKERR DMRestoreNamedGlobalVector(dummy_dm,
"TSTheta_X0",
1915 &ts_solver_vecs[0]);
1916 CHKERR DMRestoreNamedGlobalVector(dummy_dm,
"TSTheta_Xdot",
1917 &ts_solver_vecs[1]);
1920 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
1921 MOFEM_LOG(
"WORLD", Sev::verbose) <<
"Zero field " << f;
1922 CHKERR field_blas->fieldLambdaOnEntities(zero_dofs, f);
1931 auto post_proc = [&](
auto exe_test) {
1933 auto post_proc_fe = boost::make_shared<PostProcEleDomain>(
mField);
1934 auto &pip = post_proc_fe->getOpPtrVector();
1935 post_proc_fe->exeTestHook = exe_test;
1943 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1946 fe_parent->getOpPtrVector(), {H1});
1954 auto u_ptr = boost::make_shared<MatrixDouble>();
1955 auto p_ptr = boost::make_shared<VectorDouble>();
1956 auto h_ptr = boost::make_shared<VectorDouble>();
1957 auto g_ptr = boost::make_shared<VectorDouble>();
1959 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"U",
1962 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"P",
1965 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"H",
1968 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"G",
1972 post_proc_fe->getOpPtrVector().push_back(
1974 new OpPPMap(post_proc_fe->getPostProcMesh(),
1975 post_proc_fe->getMapGaussPts(),
1977 {{
"P", p_ptr}, {
"H", h_ptr}, {
"G", g_ptr}},
1989 auto dm =
simple->getDM();
1991 CHKERR post_proc_fe->writeFile(
"out_projection.h5m");
2000 if constexpr (
debug) {
2006 fe_domain_lhs.swap(pip_mng->getDomainLhsFE());
2007 fe_domain_rhs.swap(pip_mng->getDomainRhsFE());
2008 fe_bdy_lhs.swap(pip_mng->getBoundaryLhsFE());
2009 fe_bdy_rhs.swap(pip_mng->getBoundaryRhsFE());
2022 auto add_parent_field_domain = [&](
auto fe,
auto op,
auto field) {
2027 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2035 auto add_parent_field_bdy = [&](
auto fe,
auto op,
auto field) {
2040 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2053 auto dot_u_ptr = boost::make_shared<MatrixDouble>();
2054 auto u_ptr = boost::make_shared<MatrixDouble>();
2055 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
2056 auto dot_h_ptr = boost::make_shared<VectorDouble>();
2057 auto h_ptr = boost::make_shared<VectorDouble>();
2058 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
2059 auto g_ptr = boost::make_shared<VectorDouble>();
2060 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
2061 auto lambda_ptr = boost::make_shared<VectorDouble>();
2062 auto p_ptr = boost::make_shared<VectorDouble>();
2063 auto div_u_ptr = boost::make_shared<VectorDouble>();
2067 auto set_domain_general = [&](
auto fe) {
2069 auto &pip = fe->getOpPtrVector();
2077 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2080 fe_parent->getOpPtrVector(), {H1});
2086 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"U");
2106 "Coordinate system not implemented");
2109 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"H");
2115 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"G");
2120 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"P");
2125 auto set_domain_rhs = [&](
auto fe) {
2127 auto &pip = fe->getOpPtrVector();
2129 CHKERR set_domain_general(fe);
2131 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"U");
2132 pip.push_back(
new OpRhsU(
"U", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr,
2133 grad_h_ptr, g_ptr, p_ptr));
2135 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"H");
2136 pip.push_back(
new OpRhsH<false>(
"H", u_ptr, dot_h_ptr, h_ptr, grad_h_ptr,
2139 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"G");
2140 pip.push_back(
new OpRhsG<false>(
"G", h_ptr, grad_h_ptr, g_ptr));
2142 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"P");
2144 "P", div_u_ptr, [](
const double r,
const double,
const double) {
2148 "P", p_ptr, [](
const double r,
const double,
const double) {
2154 auto set_domain_lhs = [&](
auto fe) {
2156 auto &pip = fe->getOpPtrVector();
2158 CHKERR set_domain_general(fe);
2160 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"U");
2162 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"U");
2163 pip.push_back(
new OpLhsU_dU(
"U", u_ptr, grad_u_ptr, h_ptr));
2164 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"H");
2166 new OpLhsU_dH(
"U",
"H", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr, g_ptr));
2168 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"G");
2169 pip.push_back(
new OpLhsU_dG(
"U",
"G", grad_h_ptr));
2172 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"H");
2174 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"U");
2175 pip.push_back(
new OpLhsH_dU(
"H",
"U", grad_h_ptr));
2176 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"H");
2178 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"G");
2182 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"G");
2184 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"H");
2186 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"G");
2190 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"P");
2192 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"U");
2198 [](
const double r,
const double,
const double) {
2206 [](
const double r,
const double,
const double) {
2213 "Coordinate system not implemented");
2216 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"P");
2217 pip.push_back(
new OpDomainMassP(
"P",
"P", [](
double r,
double,
double) {
2225 auto get_block_name = [](
auto name) {
2226 return boost::format(
"%s(.*)") %
"WETTING_ANGLE";
2229 auto get_blocks = [&](
auto &&name) {
2231 std::regex(name.str()));
2234 auto set_boundary_rhs = [&](
auto fe) {
2236 auto &pip = fe->getOpPtrVector();
2242 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2249 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"U");
2252 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"L");
2257 pip,
mField,
"L", {},
"INFLUX",
2260 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"U");
2263 auto wetting_block = get_blocks(get_block_name(
"WETTING_ANGLE"));
2264 if (wetting_block.size()) {
2271 op_bdy_side->getOpPtrVector(), {H1});
2274 op_bdy_side->getSideFEPtr(),
"", UDO::OPSPACE,
2278 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2281 fe_parent->getOpPtrVector(), {H1});
2287 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPROW,
2289 op_bdy_side->getOpPtrVector().push_back(
2292 pip.push_back(op_bdy_side);
2295 for (
auto &b : wetting_block) {
2297 std::vector<double> attr_vec;
2298 CHKERR b->getMeshsetIdEntitiesByDimension(
2300 b->getAttributes(attr_vec);
2301 if (attr_vec.size() != 1)
2303 "Should be one attribute");
2304 MOFEM_LOG(
"FS", Sev::inform) <<
"Wetting angle: " << attr_vec.front();
2306 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"G");
2308 "G", grad_h_ptr, boost::make_shared<Range>(force_edges),
2316 auto set_boundary_lhs = [&](
auto fe) {
2318 auto &pip = fe->getOpPtrVector();
2324 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2331 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"L");
2332 CHKERR add_parent_field_bdy(fe, UDO::OPCOL,
"U");
2335 auto wetting_block = get_blocks(get_block_name(
"WETTING_ANGLE"));
2336 if (wetting_block.size()) {
2337 auto col_ind_ptr = boost::make_shared<std::vector<VectorInt>>();
2338 auto col_diff_base_ptr = boost::make_shared<std::vector<MatrixDouble>>();
2346 op_bdy_side->getOpPtrVector(), {H1});
2349 op_bdy_side->getSideFEPtr(),
"", UDO::OPSPACE,
2353 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2356 fe_parent->getOpPtrVector(), {H1});
2362 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPROW,
2364 op_bdy_side->getOpPtrVector().push_back(
2366 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPCOL,
2368 op_bdy_side->getOpPtrVector().push_back(
2372 pip.push_back(op_bdy_side);
2375 for (
auto &b : wetting_block) {
2377 std::vector<double> attr_vec;
2378 CHKERR b->getMeshsetIdEntitiesByDimension(
2380 b->getAttributes(attr_vec);
2381 if (attr_vec.size() != 1)
2383 "Should be one attribute");
2385 <<
"wetting angle edges size " << force_edges.size();
2387 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"G");
2389 "G", grad_h_ptr, col_ind_ptr, col_diff_base_ptr,
2390 boost::make_shared<Range>(force_edges), attr_vec.front()));
2399 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
2400 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
2401 CHKERR set_boundary_rhs(pip_mng->getCastBoundaryRhsFE());
2402 CHKERR set_boundary_lhs(pip_mng->getCastBoundaryLhsFE());
2426 boost::shared_ptr<PostProcEleDomainCont> post_proc,
2427 boost::shared_ptr<PostProcEleBdyCont> post_proc_edge,
2428 std::pair<boost::shared_ptr<BoundaryEle>, boost::shared_ptr<VectorDouble>>
2435 MOFEM_LOG(
"FS", Sev::verbose) <<
"Monitor";
2438 MOFEM_LOG(
"FS", Sev::verbose) <<
"Mesh pre proc";
2441 auto post_proc_begin =
2442 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
2444 auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
2452 CHKERR post_proc_end->writeFile(
2453 "out_step_" + boost::lexical_cast<std::string>(
ts_step) +
".h5m");
2454 MOFEM_LOG(
"FS", Sev::verbose) <<
"Mesh pre proc done";
2460 MPI_Allreduce(MPI_IN_PLACE, &(*
liftVec)[0],
SPACE_DIM, MPI_DOUBLE, MPI_SUM,
2464 <<
" lift vec x: " << (*liftVec)[0] <<
" y: " << (*liftVec)[1];
2490 auto setup_subdm = [&](
auto dm) {
2494 auto dm =
simple->getDM();
2495 auto level_ents_ptr = boost::make_shared<Range>();
2496 CHKERR bit_mng->getEntitiesByRefLevel(
2499 CHKERR DMSetType(subdm,
"DMMOFEM");
2505 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
2518 auto get_fe_post_proc = [&](
auto post_proc_mesh) {
2519 auto add_parent_field_domain = [&](
auto fe,
auto op,
auto field) {
2524 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2533 boost::make_shared<PostProcEleDomainCont>(
mField, post_proc_mesh);
2534 post_proc_fe->exeTestHook = [](
FEMethod *fe_ptr) {
2539 auto u_ptr = boost::make_shared<MatrixDouble>();
2540 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
2541 auto h_ptr = boost::make_shared<VectorDouble>();
2542 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
2543 auto p_ptr = boost::make_shared<VectorDouble>();
2544 auto g_ptr = boost::make_shared<VectorDouble>();
2545 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
2548 post_proc_fe->getOpPtrVector(), {H1});
2554 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2557 fe_parent->getOpPtrVector(), {H1});
2563 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"U");
2564 post_proc_fe->getOpPtrVector().push_back(
2566 post_proc_fe->getOpPtrVector().push_back(
2570 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"H");
2571 post_proc_fe->getOpPtrVector().push_back(
2573 post_proc_fe->getOpPtrVector().push_back(
2576 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"P");
2577 post_proc_fe->getOpPtrVector().push_back(
2580 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"G");
2581 post_proc_fe->getOpPtrVector().push_back(
2583 post_proc_fe->getOpPtrVector().push_back(
2588 post_proc_fe->getOpPtrVector().push_back(
2591 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
2593 {{
"H", h_ptr}, {
"P", p_ptr}, {
"G", g_ptr}},
2595 {{
"U", u_ptr}, {
"H_GRAD", grad_h_ptr}, {
"G_GRAD", grad_g_ptr}},
2597 {{
"GRAD_U", grad_u_ptr}},
2605 return post_proc_fe;
2608 auto get_bdy_post_proc_fe = [&](
auto post_proc_mesh) {
2609 auto add_parent_field_bdy = [&](
auto fe,
auto op,
auto field) {
2610 return setParentDofs(
2614 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2623 boost::make_shared<PostProcEleBdyCont>(mField, post_proc_mesh);
2624 post_proc_fe->exeTestHook = [](
FEMethod *fe_ptr) {
2633 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2642 OpGetNormal(boost::shared_ptr<VectorDouble> l_ptr,
2643 boost::shared_ptr<MatrixDouble> n_ptr)
2651 auto t_n_fe = getFTensor1NormalsAtGaussPts();
2652 ptrNormal->resize(
SPACE_DIM, getGaussPts().size2(),
false);
2653 auto t_n = getFTensor1FromMat<SPACE_DIM>(*ptrNormal);
2654 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
2655 t_n(
i) = t_n_fe(
i) * t_l / std::sqrt(t_n_fe(
i) * t_n_fe(
i));
2664 boost::shared_ptr<VectorDouble> ptrL;
2665 boost::shared_ptr<MatrixDouble> ptrNormal;
2668 auto u_ptr = boost::make_shared<MatrixDouble>();
2669 auto p_ptr = boost::make_shared<VectorDouble>();
2670 auto lambda_ptr = boost::make_shared<VectorDouble>();
2671 auto normal_l_ptr = boost::make_shared<MatrixDouble>();
2673 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW,
"U");
2674 post_proc_fe->getOpPtrVector().push_back(
2677 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW,
"L");
2678 post_proc_fe->getOpPtrVector().push_back(
2680 post_proc_fe->getOpPtrVector().push_back(
2681 new OpGetNormal(lambda_ptr, normal_l_ptr));
2683 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW,
"P");
2684 post_proc_fe->getOpPtrVector().push_back(
2688 op_ptr->doWorkRhsHook = [&](
DataOperator *base_op_ptr,
int side,
2697 post_proc_fe->getOpPtrVector().push_back(
2699 new OpPPMap(post_proc_fe->getPostProcMesh(),
2700 post_proc_fe->getMapGaussPts(),
2714 return post_proc_fe;
2717 auto get_lift_fe = [&]() {
2718 auto add_parent_field_bdy = [&](
auto fe,
auto op,
auto field) {
2719 return setParentDofs(
2723 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2731 auto fe = boost::make_shared<BoundaryEle>(mField);
2732 fe->exeTestHook = [](
FEMethod *fe_ptr) {
2737 auto lift_ptr = boost::make_shared<VectorDouble>();
2738 auto p_ptr = boost::make_shared<VectorDouble>();
2739 auto ents_ptr = boost::make_shared<Range>();
2745 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2752 std::vector<const CubitMeshSets *> vec_ptr;
2754 std::regex(
"LIFT"), vec_ptr);
2755 for (
auto m_ptr : vec_ptr) {
2760 ents_ptr->merge(ents);
2763 MOFEM_LOG(
"FS", Sev::noisy) <<
"Lift ents " << (*ents_ptr);
2765 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"P");
2766 fe->getOpPtrVector().push_back(
2768 fe->getOpPtrVector().push_back(
2771 return std::make_pair(fe, lift_ptr);
2774 auto set_post_proc_monitor = [&](
auto ts) {
2778 boost::shared_ptr<FEMethod> null_fe;
2779 auto post_proc_mesh = boost::make_shared<moab::Core>();
2780 auto monitor_ptr = boost::make_shared<Monitor>(
2782 get_fe_post_proc(post_proc_mesh), get_bdy_post_proc_fe(post_proc_mesh),
2785 null_fe, monitor_ptr);
2789 auto dm =
simple->getDM();
2793 auto ts_pre_post_proc = boost::make_shared<TSPrePostProc>();
2798 ptr->fsRawPtr =
this;
2799 ptr->solverSubDM = create_solver_dm(
simple->getDM());
2803 CHKERR VecAssemblyBegin(ptr->globSol);
2804 CHKERR VecAssemblyEnd(ptr->globSol);
2806 auto sub_ts = pip_mng->createTSIM(ptr->solverSubDM);
2808 CHKERR set_post_proc_monitor(sub_ts);
2811 CHKERR TSSetFromOptions(ts);
2815 auto print_fields_in_section = [&]() {
2818 simple->getProblemName());
2819 PetscInt num_fields;
2820 CHKERR PetscSectionGetNumFields(section, &num_fields);
2821 for (
int f = 0;
f < num_fields; ++
f) {
2825 <<
"Field " <<
f <<
" " << std::string(
field_name);
2830 CHKERR print_fields_in_section();
2832 CHKERR TSSolve(ts, ptr->globSol);
2857#ifdef PYTHON_INIT_SURFACE
2862 const char param_file[] =
"param_file.petsc";
2866 auto core_log = logging::core::get();
2874 DMType dm_name =
"DMMOFEM";
2879 moab::Core mb_instance;
2880 moab::Interface &moab = mb_instance;
2897#ifdef PYTHON_INIT_SURFACE
2898 if (Py_FinalizeEx() < 0) {
2914 "can not get vertices on bit 0");
2919 Range plus_range, minus_range;
2920 std::vector<EntityHandle> plus, minus;
2923 for (
auto p = vertices.pair_begin(); p != vertices.pair_end(); ++p) {
2925 const auto f = p->first;
2926 const auto s = p->second;
2931 auto it = dofs_mi.lower_bound(lo_uid);
2932 const auto hi_it = dofs_mi.upper_bound(hi_uid);
2936 plus.reserve(std::distance(it, hi_it));
2937 minus.reserve(std::distance(it, hi_it));
2939 for (; it != hi_it; ++it) {
2940 const auto v = (*it)->getFieldData();
2942 plus.push_back((*it)->getEnt());
2944 minus.push_back((*it)->getEnt());
2947 plus_range.insert_list(plus.begin(), plus.end());
2948 minus_range.insert_list(minus.begin(), minus.end());
2953 <<
"Plus range " << plus_range << endl;
2955 <<
"Minus range " << minus_range << endl;
2958 auto get_elems = [&](
auto &ents,
auto bit,
auto mask) {
2961 moab::Interface::UNION),
2962 "can not get adjacencies");
2964 "can not filter elements with bit 0");
2968 CHKERR comm_mng->synchroniseEntities(plus_range);
2969 CHKERR comm_mng->synchroniseEntities(minus_range);
2972 ele_plus[0] = get_elems(plus_range,
bit(0),
BitRefLevel().set());
2973 ele_minus[0] = get_elems(minus_range,
bit(0),
BitRefLevel().set());
2974 auto common = intersect(ele_plus[0], ele_minus[0]);
2975 ele_plus[0] = subtract(ele_plus[0], common);
2976 ele_minus[0] = subtract(ele_minus[0], common);
2978 auto get_children = [&](
auto &p,
auto &
c) {
2980 CHKERR bit_mng->updateRangeByChildren(p,
c);
2992 auto get_level = [&](
auto &p,
auto &
m,
auto z,
auto bit,
auto mask) {
2995 bit_mng->getEntitiesByDimAndRefLevel(
bit, mask,
SPACE_DIM,
l),
2996 "can not get vertices on bit");
2999 for (
auto f = 0; f != z; ++f) {
3003 l = get_elems(conn,
bit, mask);
3008 std::vector<Range> vec_levels(
nb_levels);
3010 vec_levels[
l] = get_level(ele_plus[
l], ele_minus[
l], 2 * overlap,
bit(
l),
3014 if constexpr (
debug) {
3017 std::string name = (boost::format(
"out_r%d.h5m") %
l).str();
3034 start_mask[s] =
true;
3040 Range prev_level_skin;
3044 CHKERR bit_mng->lambdaBitRefLevel(
3051 auto set_levels = [&](
auto &&
3061 auto get_adj = [&](
auto ents) {
3067 ents.subset_by_dimension(
SPACE_DIM), d,
false, conn,
3068 moab::Interface::UNION),
3082 CHKERR bit_mng->updateRangeByParent(vec_levels[
l], parents);
3085 level_prev = subtract(level_prev, parents);
3087 level_prev.merge(vec_levels[
l]);
3095 CHKERR bit_mng->getEntitiesByDimAndRefLevel(
3097 level = get_adj(level);
3106 auto set_skins = [&]() {
3110 ParallelComm *pcomm =
3119 "can't get bit level");
3126 auto get_level_skin = [&]() {
3132 CHKERR pcomm->filter_pstatus(skin_level_mesh,
3133 PSTATUS_SHARED | PSTATUS_MULTISHARED,
3134 PSTATUS_NOT, -1,
nullptr);
3137 skin_level = subtract(skin_level,
3141 Range skin_level_verts;
3144 skin_level.merge(skin_level_verts);
3147 bit_prev.set(
l - 1);
3151 skin.merge(subtract(skin_level, level_prev));
3157 auto resolve_shared = [&](
auto &&skin) {
3158 Range tmp_skin = skin;
3160 map<int, Range> map_procs;
3162 tmp_skin, &map_procs);
3164 Range from_other_procs;
3165 for (
auto &
m : map_procs) {
3167 from_other_procs.merge(
m.second);
3171 auto common = intersect(
3172 skin, from_other_procs);
3173 skin.merge(from_other_procs);
3177 if (!common.empty()) {
3179 skin = subtract(skin, common);
3186 auto get_parent_level_skin = [&](
auto skin) {
3188 CHKERR bit_mng->updateRangeByParent(
3189 skin.subset_by_dimension(
SPACE_DIM - 1), skin_parents);
3190 Range skin_parent_verts;
3193 skin_parents.merge(skin_parent_verts);
3196 return skin_parents;
3199 auto child_skin = resolve_shared(get_level_skin());
3200 auto parent_skin = get_parent_level_skin(child_skin);
3202 child_skin = subtract(child_skin, parent_skin);
3210 auto set_current = [&]() {
3219 last_level = subtract(last_level, skin_child);
3231 if constexpr (
debug) {
3234 std::string name = (boost::format(
"out_level%d.h5m") %
l).str();
3237 "PARALLEL=WRITE_PART");
3241 "MOAB",
"PARALLEL=WRITE_PART");
3244 "MOAB",
"PARALLEL=WRITE_PART");
3248 "MOAB",
"PARALLEL=WRITE_PART");
3251 "MOAB",
"PARALLEL=WRITE_PART");
3259 boost::shared_ptr<FEMethod> fe_top, std::string
field_name,
3260 ForcesAndSourcesCore::UserDataOperator::OpType op,
3262 boost::function<boost::shared_ptr<ForcesAndSourcesCore>()> get_elem,
3269 boost::function<void(boost::shared_ptr<ForcesAndSourcesCore>,
int)>
3271 [&](boost::shared_ptr<ForcesAndSourcesCore> parent_fe_pt,
int level) {
3276 auto fe_ptr_current = get_elem();
3280 boost::dynamic_pointer_cast<ForcesAndSourcesCore>(
3285 if (op == ForcesAndSourcesCore::UserDataOperator::OPSPACE) {
3288 parent_fe_pt->getOpPtrVector().push_back(
3292 H1, op, fe_ptr_current,
3303 parent_fe_pt->getOpPtrVector().push_back(
3318 add_parent_level(boost::dynamic_pointer_cast<ForcesAndSourcesCore>(fe_top),
3333 auto cut_off_dofs = [&]() {
3336 auto &m_field = ptr->fsRawPtr->mField;
3338 Range current_verts;
3343 auto cut_off_verts = [&](boost::shared_ptr<FieldEntity> ent_ptr) {
3345 for (
auto &
h : ent_ptr->getEntFieldData()) {
3351 auto field_blas = m_field.getInterface<
FieldBlas>();
3361 MOFEM_LOG(
"FS", Sev::inform) <<
"Run step pre proc";
3363 auto &m_field = ptr->fsRawPtr->mField;
3367 auto get_norm = [&](
auto x) {
3369 CHKERR VecNorm(x, NORM_2, &nrm);
3374 auto refine_problem = [&]() {
3376 MOFEM_LOG(
"FS", Sev::inform) <<
"Refine problem";
3378 CHKERR ptr->fsRawPtr->projectData();
3384 auto set_jacobian_operators = [&]() {
3387 CHKERR KSPReset(ptr->subKSP);
3392 auto set_solution = [&]() {
3394 MOFEM_LOG(
"FS", Sev::inform) <<
"Set solution";
3396 PetscObjectState state;
3407 INSERT_VALUES, SCATTER_FORWARD);
3410 <<
"Set solution, vector norm " << get_norm(ptr->globSol);
3415 PetscObjectTypeCompare((PetscObject)ts, TSTHETA, &is_theta);
3419 CHKERR set_jacobian_operators();
3424 "Sorry, only TSTheta handling is implemented");
3429 PetscBarrier((PetscObject)ts);
3441 auto &m_field = ptr->fsRawPtr->mField;
3453 auto sub_u = ptr->getSubVector();
3456 auto scatter = ptr->getScatter(sub_u, u,
R);
3458 auto apply_scatter_and_update = [&](
auto x,
auto sub_x) {
3460 CHKERR VecScatterBegin(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3461 CHKERR VecScatterEnd(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3462 CHKERR VecGhostUpdateBegin(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3463 CHKERR VecGhostUpdateEnd(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3467 CHKERR apply_scatter_and_update(u, sub_u);
3468 CHKERR apply_scatter_and_update(u_t, sub_u_t);
3471 CHKERR VecScatterBegin(scatter, sub_f, f, INSERT_VALUES, SCATTER_FORWARD);
3472 CHKERR VecScatterEnd(scatter, sub_f, f, INSERT_VALUES, SCATTER_FORWARD);
3478 PetscReal
a, Mat
A, Mat
B,
3482 auto sub_u = ptr->getSubVector();
3484 auto scatter = ptr->getScatter(sub_u, u,
R);
3486 auto apply_scatter_and_update = [&](
auto x,
auto sub_x) {
3488 CHKERR VecScatterBegin(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3489 CHKERR VecScatterEnd(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3490 CHKERR VecGhostUpdateBegin(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3491 CHKERR VecGhostUpdateEnd(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3495 CHKERR apply_scatter_and_update(u, sub_u);
3496 CHKERR apply_scatter_and_update(u_t, sub_u_t);
3499 ptr->tsCtxPtr.get());
3508 auto get_norm = [&](
auto x) {
3510 CHKERR VecNorm(x, NORM_2, &nrm);
3514 auto sub_u = ptr->getSubVector();
3515 auto scatter = ptr->getScatter(sub_u, u,
R);
3516 CHKERR VecScatterBegin(scatter, u, sub_u, INSERT_VALUES, SCATTER_REVERSE);
3517 CHKERR VecScatterEnd(scatter, u, sub_u, INSERT_VALUES, SCATTER_REVERSE);
3518 CHKERR VecGhostUpdateBegin(sub_u, INSERT_VALUES, SCATTER_FORWARD);
3519 CHKERR VecGhostUpdateEnd(sub_u, INSERT_VALUES, SCATTER_FORWARD);
3522 <<
"u norm " << get_norm(u) <<
" u sub nom " << get_norm(sub_u);
3532 MOFEM_LOG(
"FS", Sev::verbose) <<
"SetUP sub PC";
3533 ptr->subKSP =
createKSP(ptr->fsRawPtr->mField.get_comm());
3534 CHKERR KSPSetFromOptions(ptr->subKSP);
3542 auto sub_x = ptr->getSubVector();
3544 auto scatter = ptr->getScatter(sub_x, pc_x,
R);
3545 CHKERR VecScatterBegin(scatter, pc_f, sub_f, INSERT_VALUES,
3547 CHKERR VecScatterEnd(scatter, pc_f, sub_f, INSERT_VALUES, SCATTER_REVERSE);
3548 CHKERR KSPSetOperators(ptr->subKSP, ptr->subB, ptr->subB);
3549 MOFEM_LOG(
"FS", Sev::verbose) <<
"PCShell solve";
3550 CHKERR KSPSolve(ptr->subKSP, sub_f, sub_x);
3551 MOFEM_LOG(
"FS", Sev::verbose) <<
"PCShell solve <- done";
3552 CHKERR VecScatterBegin(scatter, sub_x, pc_x, INSERT_VALUES,
3554 CHKERR VecScatterEnd(scatter, sub_x, pc_x, INSERT_VALUES, SCATTER_FORWARD);
3561 auto prb_ptr = ptr->fsRawPtr->mField.get_problem(
"SUB_SOLVER");
3562 if (
auto sub_data = prb_ptr->getSubData()) {
3563 auto is = sub_data->getSmartColIs();
3598 CHKERR TSGetSNES(ts, &snes);
3601 auto set_section_monitor = [&](
auto snes) {
3603 CHKERR SNESMonitorSet(snes,
3606 (
void *)(snes_ctx_ptr.get()),
nullptr);
3610 CHKERR set_section_monitor(snes);
3612 auto ksp =
createKSP(m_field.get_comm());
3613 CHKERR KSPSetType(ksp, KSPPREONLY);
3615 CHKERR PCSetType(sub_pc, PCSHELL);
3618 CHKERR KSPSetPC(ksp, sub_pc);
3619 CHKERR SNESSetKSP(snes, ksp);
#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)
ElementsAndOps< SPACE_DIM >::DomainParentEle DomainParentEle
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
BoundaryEle::UserDataOperator BoundaryEleOp
Kronecker Delta class symmetric.
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define BITREFLEVEL_SIZE
max number of refinements
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#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 ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
auto cylindrical
Coordinate system scaling factor.
auto marker
Create marker bit reference level
auto bubble_device
Bubble device initialization.
ElementsAndOps< SPACE_DIM >::DomainParentEle DomainParentEle
auto save_range
Save range of entities to file.
auto get_M2
Degenerate mobility function M₂
constexpr int U_FIELD_DIM
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
static boost::weak_ptr< TSPrePostProc > tsPrePostProc
FTensor::Index< 'j', SPACE_DIM > j
auto get_M2_dh
Derivative of degenerate mobility M₂
auto kernel_eye
Eye-shaped interface initialization
double lambda
surface tension
FTensor::Index< 'k', SPACE_DIM > k
OpDomainMassH OpDomainMassG
auto get_M3_dh
Derivative of non-linear mobility M₃
FTensor::Index< 'i', SPACE_DIM > i
auto init_h
Initialisation function.
auto d_phase_function_h
Derivative of phase function with respect to h.
auto get_fe_bit
Get bit reference level from finite element.
FTensor::Index< 'l', SPACE_DIM > l
auto kernel_oscillation
Oscillating interface initialization.
auto get_M0
Constant mobility function M₀
auto get_f_dh
Derivative of double-well potential.
auto get_M3
Non-linear mobility function M₃
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesVector< BASE_DIM, SPACE_DIM, 1 > OpDomainAssembleVector
OpDomainMassH OpDomainMassP
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, 1 > OpDomainMassH
auto wetting_angle
Wetting angle function (placeholder)
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpDomainAssembleScalar
auto get_dofs_ents_all
Get all entities with DOFs in the problem - used for debugging.
auto phase_function
Phase-dependent material property interpolation.
constexpr IntegrationType I
auto get_skin_projection_bit
auto get_global_size
Get global size across all processors.
auto get_dofs_ents_by_field_name
Get entities of DOFs by field name - used for debugging.
auto get_D
Create deviatoric stress tensor.
auto wetting_angle_sub_stepping
Wetting angle sub-stepping for gradual application.
auto get_M0_dh
Derivative of constant mobility.
auto my_max
Maximum function with smooth transition.
FormsIntegrators< BoundaryEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, 1 > OpBoundaryMassL
auto cut_off
Phase field cutoff function.
int order
approximation order
SideEle::UserDataOperator SideOp
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMixScalarTimesDiv< SPACE_DIM, COORD_TYPE > OpMixScalarTimesDiv
FormsIntegrators< BoundaryEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpBoundaryAssembleScalar
auto get_current_bit
dofs bit used to do calculations
BoundaryEle::UserDataOperator BoundaryEleOp
auto bit
Create bit reference level.
auto capillary_tube
Capillary tube initialization.
auto get_f
Double-well potential function.
auto d_cut_off
Derivative of cutoff function.
auto my_min
Minimum function with smooth transition
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, U_FIELD_DIM > OpDomainMassU
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
auto createDMVector(DM dm)
Get smart vector from DM.
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
auto createDMMatrix(DM dm)
Get smart matrix from DM.
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
virtual const DofEntity_multiIndex * get_dofs() const =0
Get the dofs object.
MoFEMErrorCode getEntitiesByDimAndRefLevel(const BitRefLevel bit, const BitRefLevel mask, const int dim, const EntityHandle meshset, int verb=0) const
add all ents from ref level given by bit to meshset
MoFEMErrorCode getEntitiesByTypeAndRefLevel(const BitRefLevel bit, const BitRefLevel mask, const EntityType type, const EntityHandle meshset, int verb=0) const
add all ents from ref level given by bit to meshset
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 getMeshset(const int ms_id, const unsigned int cubit_bc_type, EntityHandle &meshset) const
get meshset from CUBIT Id and CUBIT type
MoFEMErrorCode buildSubProblem(const std::string out_name, const std::vector< std::string > &fields_row, const std::vector< std::string > &fields_col, const std::string main_problem, const bool square_matrix=true, const map< std::string, boost::shared_ptr< Range > > *entityMapRow=nullptr, const map< std::string, boost::shared_ptr< Range > > *entityMapCol=nullptr, int verb=VERBOSE)
build sub problem
auto test_bit_child
lambda function used to select elements on which finite element pipelines are executed.
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
implementation of Data Operators for Forces and Sources
auto createKSP(MPI_Comm comm)
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.
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobian in TS solver.
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
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.
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
PetscErrorCode TsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t, Vec F, void *ctx)
Set IFunction for TS solver.
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto createVecScatter(Vec x, IS ix, Vec y, IS iy)
Create a Vec Scatter object.
static PetscErrorCode solve(Mat mat, Vec x, Vec y)
auto createTS(MPI_Comm comm)
auto createPC(MPI_Comm comm)
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
PetscObject getPetscObject(T obj)
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
auto getProblemPtr(DM dm)
get problem pointer from DM
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
ElementsAndOps< SPACE_DIM >::SideEle SideEle
constexpr double t
plate stiffness
constexpr auto field_name
OpBaseImpl< PETSC, EdgeEleOp > OpBase
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
FTensor::Index< 'm', 3 > m
MoFEM::Interface & mField
Reference to MoFEM interface.
Range findParentsToRefine(Range ents, BitRefLevel level, BitRefLevel mask)
Find parent entities that need refinement.
MoFEMErrorCode setupProblem()
Setup problem fields and parameters.
MoFEMErrorCode runProblem()
Main function to run the complete free surface simulation.
FreeSurface(MoFEM::Interface &m_field)
Constructor.
MoFEMErrorCode projectData()
Project solution data between mesh levels.
MoFEMErrorCode boundaryCondition()
Apply boundary conditions and initialize fields.
MoFEMErrorCode readMesh()
Read mesh from input file.
MoFEMErrorCode setParentDofs(boost::shared_ptr< FEMethod > fe_top, std::string field_name, ForcesAndSourcesCore::UserDataOperator::OpType op, BitRefLevel child_ent_bit, boost::function< boost::shared_ptr< ForcesAndSourcesCore >()> get_elem, int verbosity, LogManager::SeverityLevel sev)
Create hierarchy of elements run on parent levels.
std::vector< Range > findEntitiesCrossedByPhaseInterface(size_t overlap)
Find entities on refinement levels.
MoFEMErrorCode assembleSystem()
Assemble system operators and matrices.
MoFEM::Interface & mField
MoFEMErrorCode refineMesh(size_t overlap)
Perform adaptive mesh refinement.
MoFEMErrorCode makeRefProblem()
Create refined problem for mesh adaptation.
MoFEMErrorCode solveSystem()
Solve the time-dependent free surface problem.
Add operators pushing bases from local to physical configuration.
boost::weak_ptr< CacheTuple > getCacheWeakPtr() const
Get the cache weak pointer object.
Boundary condition manager for finite element problem setup.
virtual int get_comm_size() const =0
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
base operator to do operations at Gauss Pt. level
Deprecated interface functions.
static UId getLoFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
static UId getHiFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
Data on single entity (This is passed as argument to DataOperator::doWork)
Structure for user loop methods on finite elements.
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
Shared pointer to finite element database structure.
MoFEMErrorCode fieldLambdaOnEntities(OneFieldFunctionOnEntities lambda, const std::string field_name, Range *ents_ptr=nullptr)
field lambda
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
Section manager is used to create indexes and sections.
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.
Interface for managing meshsets containing materials and boundary conditions.
Operator to project base functions from parent entity to child.
Calculate divergence of vector field at integration points.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Specialization for double precision scalar field values calculation.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Approximate field values for given petsc vector.
Specialization for double precision vector field values calculation.
Element used to execute operators on side of the element.
Post post-proc data at points from hash maps.
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
Operator to execute finite element instance on parent element. This operator is typically used to pro...
Template struct for dimension-specific finite element types.
PipelineManager interface.
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
MoFEMErrorCode getDM(DM *dm)
Get DM.
bool & getParentAdjacencies()
Get the addParentAdjacencies flag.
intrusive_ptr for managing petsc objects
PetscReal ts_t
Current time value.
PetscInt ts_step
Current time step number.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
[Push operators to pipeline]
boost::shared_ptr< moab::Core > postProcMesh
MoFEMErrorCode postProcess()
Post-processing function executed at loop completion.
Monitor(SmartPetscObj< DM > dm, boost::shared_ptr< moab::Core > post_proc_mesh, boost::shared_ptr< PostProcEleDomainCont > post_proc, boost::shared_ptr< PostProcEleBdyCont > post_proc_edge, std::pair< boost::shared_ptr< BoundaryEle >, boost::shared_ptr< VectorDouble > > p)
boost::shared_ptr< VectorDouble > liftVec
boost::shared_ptr< BoundaryEle > liftFE
boost::shared_ptr< PostProcEle > postProc
boost::shared_ptr< PostProcEleBdyCont > postProcEdge
boost::shared_ptr< PostProcEleDomainCont > postProc
Set of functions called by PETSc solver used to refine and update mesh.
boost::shared_ptr< TsCtx > tsCtxPtr
virtual ~TSPrePostProc()=default
SmartPetscObj< Vec > globRes
SmartPetscObj< Mat > subB
SmartPetscObj< Vec > globSol
SmartPetscObj< Vec > getSubVector()
Create sub-problem vector.
static MoFEMErrorCode tsMonitor(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Monitor solution during time stepping.
static MoFEMErrorCode pcApply(PC pc, Vec pc_f, Vec pc_x)
Apply preconditioner.
SmartPetscObj< VecScatter > getScatter(Vec x, Vec y, enum FR fr)
Get scatter context for vector operations.
SmartPetscObj< DM > solverSubDM
static MoFEMErrorCode tsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set implicit Jacobian for time stepping.
boost::shared_ptr< SnesCtx > snesCtxPtr
static MoFEMErrorCode tsPreProc(TS ts)
Pre process time step.
static MoFEMErrorCode tsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t, Vec f, void *ctx)
Set implicit function for time stepping.
static MoFEMErrorCode pcSetup(PC pc)
Setup preconditioner.
static MoFEMErrorCode tsPostProc(TS ts)
Post process time step.
static MoFEMErrorCode tsPreStage(TS ts)
Pre-stage processing for time stepping.
MoFEMErrorCode tsSetUp(TS ts)
Used to setup TS solver.
SmartPetscObj< KSP > subKSP
constexpr CoordinateTypes COORD_TYPE