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;
35 auto main_module = bp::import(
"__main__");
36 mainNamespace = main_module.attr(
"__dict__");
37 bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
39 surfaceFun = mainNamespace[
"surface"];
41 }
catch (bp::error_already_set
const &) {
51 double x,
double y,
double z,
double eta,
double &s
58 s = bp::extract<double>(surfaceFun(x, y, z,
eta));
60 }
catch (bp::error_already_set
const &) {
69 bp::object mainNamespace;
70 bp::object surfaceFun;
73static boost::weak_ptr<SurfacePython> surfacePythonWeakPtr;
85 IntegrationType::GAUSS;
97using SideOp = SideEle::UserDataOperator;
125template <CoordinateTypes COORD_TYPE>
178double tol = std::numeric_limits<float>::epsilon();
197 constexpr int sub_stepping = 16;
198 return std::min(1.,
static_cast<double>(ts_step) / sub_stepping);
203auto my_max = [](
const double x) {
return (x - 1 + std::abs(x + 1)) / 2; };
204auto my_min = [](
const double x) {
return (x + 1 - std::abs(x - 1)) / 2; };
207 if (
h >= -1 &&
h < 1)
214 return diff *
h + ave;
223auto get_f = [](
const double h) {
return 4 *
W *
h * (
h *
h - 1); };
230 return md * (1 -
h *
h);
236 const double h2 =
h *
h;
237 const double h3 = h2 *
h;
239 return md * (2 * h3 - 3 * h2 + 1);
241 return md * (-2 * h3 - 3 * h2 + 1);
246 return md * (6 *
h * (
h - 1));
248 return md * (-6 *
h * (
h + 1));
264 constexpr double R = 0.0125;
265 constexpr double A =
R * 0.2;
266 const double theta = atan2(r, y);
267 const double w =
R + A * cos(
n * theta);
268 const double d = std::sqrt(r * r + y * y);
269 return tanh((w - d) / (
eta * std::sqrt(2)));
273 constexpr double y0 = 0.5;
274 constexpr double R = 0.4;
276 const double d = std::sqrt(r * r + y * y);
277 return tanh((
R - d) / (
eta * std::sqrt(2)));
281 constexpr double water_height = 0.;
282 return tanh((water_height - y) / (
eta * std::sqrt(2)));
287 return -tanh((-0.039 - x) / (
eta * std::sqrt(2)));
297auto init_h = [](
double r,
double y,
double theta) {
298#ifdef PYTHON_INIT_SURFACE
300 if (
auto ptr = surfacePythonWeakPtr.lock()) {
302 "error eval python");
317 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel();
322 MPI_Allreduce(&l_size, &g_size, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
326auto save_range = [](moab::Interface &moab,
const std::string name,
331 CHKERR moab.add_entities(*out_meshset, r);
332 CHKERR moab.write_file(name.c_str(),
"MOAB",
"PARALLEL=WRITE_PART",
333 out_meshset->get_ptr(), 1);
344 std::vector<EntityHandle> ents_vec;
350 auto dofs = prb_ptr->numeredRowDofsPtr;
353 ents_vec.reserve(std::distance(lo_it, hi_it));
355 for (; lo_it != hi_it; ++lo_it) {
356 ents_vec.push_back((*lo_it)->getEnt());
359 std::sort(ents_vec.begin(), ents_vec.end());
360 auto it = std::unique(ents_vec.begin(), ents_vec.end());
362 r.insert_list(ents_vec.begin(), it);
372 std::vector<EntityHandle> ents_vec;
377 auto dofs = prb_ptr->numeredRowDofsPtr;
378 ents_vec.reserve(dofs->size());
380 for (
auto d : *dofs) {
381 ents_vec.push_back(d->getEnt());
384 std::sort(ents_vec.begin(), ents_vec.end());
385 auto it = std::unique(ents_vec.begin(), ents_vec.end());
387 r.insert_list(ents_vec.begin(), it);
452 PetscReal
a, Mat A, Mat
B,
463 boost::shared_ptr<SnesCtx>
465 boost::shared_ptr<TsCtx>
514 boost::shared_ptr<FEMethod> fe_top, std::string
field_name,
515 ForcesAndSourcesCore::UserDataOperator::OpType op,
517 boost::function<boost::shared_ptr<ForcesAndSourcesCore>()> get_elem,
538 MOFEM_LOG(
"FS", Sev::inform) <<
"Read mesh for problem";
561 const char *coord_type_names[] = {
"cartesian",
"polar",
"cylindrical",
566 MOFEM_LOG(
"FS", Sev::inform) <<
"Approximation order = " <<
order;
568 <<
"Number of refinement levels nb_levels = " <<
nb_levels;
575 "-rho_m", &
rho_m, PETSC_NULLPTR);
577 &
mu_m, PETSC_NULLPTR);
579 &
rho_p, PETSC_NULLPTR);
581 &
mu_p, PETSC_NULLPTR);
585 "Height of the well in energy functional",
"-W",
599 MOFEM_LOG(
"FS", Sev::inform) <<
"Acceleration a0 = " <<
a0;
600 MOFEM_LOG(
"FS", Sev::inform) <<
"\"Minus\" phase density rho_m = " <<
rho_m;
601 MOFEM_LOG(
"FS", Sev::inform) <<
"\"Minus\" phase viscosity mu_m = " <<
mu_m;
602 MOFEM_LOG(
"FS", Sev::inform) <<
"\"Plus\" phase density rho_p = " <<
rho_p;
603 MOFEM_LOG(
"FS", Sev::inform) <<
"\"Plus\" phase viscosity mu_p = " <<
mu_p;
606 <<
"Height of the well in energy functional W = " <<
W;
607 MOFEM_LOG(
"FS", Sev::inform) <<
"Characteristic mesh size h = " <<
h;
608 MOFEM_LOG(
"FS", Sev::inform) <<
"Mobility md = " <<
md;
612 MOFEM_LOG(
"FS", Sev::inform) <<
"Average viscosity mu_ave = " <<
mu_ave;
613 MOFEM_LOG(
"FS", Sev::inform) <<
"Difference viscosity mu_diff = " <<
mu_diff;
643 auto set_problem_bit = [&]() {
666#ifdef PYTHON_INIT_SURFACE
667 auto get_py_surface_init = []() {
668 auto py_surf_init = boost::make_shared<SurfacePython>();
669 CHKERR py_surf_init->surfaceInit(
"surface.py");
670 surfacePythonWeakPtr = py_surf_init;
673 auto py_surf_init = get_py_surface_init();
680 auto dm =
simple->getDM();
682 using UDO = ForcesAndSourcesCore::UserDataOperator;
684 auto reset_bits = [&]() {
688 start_mask[s] =
true;
690 CHKERR bit_mng->lambdaBitRefLevel(
699 auto add_parent_field = [&](
auto fe,
auto op,
auto field) {
704 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
712 auto h_ptr = boost::make_shared<VectorDouble>();
713 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
714 auto g_ptr = boost::make_shared<VectorDouble>();
715 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
717 auto set_generic = [&](
auto fe) {
719 auto &pip = fe->getOpPtrVector();
727 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
730 fe_parent->getOpPtrVector(), {H1});
736 CHKERR add_parent_field(fe, UDO::OPROW,
"H");
741 CHKERR add_parent_field(fe, UDO::OPROW,
"G");
749 auto post_proc = [&](
auto exe_test) {
751 auto post_proc_fe = boost::make_shared<PostProcEleDomain>(
mField);
752 post_proc_fe->exeTestHook = exe_test;
754 CHKERR set_generic(post_proc_fe);
758 post_proc_fe->getOpPtrVector().push_back(
760 new OpPPMap(post_proc_fe->getPostProcMesh(),
761 post_proc_fe->getMapGaussPts(),
763 {{
"H", h_ptr}, {
"G", g_ptr}},
765 {{
"GRAD_H", grad_h_ptr}, {
"GRAD_G", grad_g_ptr}},
776 CHKERR post_proc_fe->writeFile(
"out_init.h5m");
781 auto solve_init = [&](
auto exe_test) {
784 pip_mng->getOpDomainRhsPipeline().clear();
785 pip_mng->getOpDomainLhsPipeline().clear();
787 auto set_domain_rhs = [&](
auto fe) {
790 auto &pip = fe->getOpPtrVector();
792 CHKERR add_parent_field(fe, UDO::OPROW,
"H");
793 pip.push_back(
new OpRhsH<true>(
"H",
nullptr,
nullptr, h_ptr, grad_h_ptr,
795 CHKERR add_parent_field(fe, UDO::OPROW,
"G");
796 pip.push_back(
new OpRhsG<true>(
"G", h_ptr, grad_h_ptr, g_ptr));
800 auto set_domain_lhs = [&](
auto fe) {
804 auto &pip = fe->getOpPtrVector();
806 CHKERR add_parent_field(fe, UDO::OPROW,
"H");
807 CHKERR add_parent_field(fe, UDO::OPCOL,
"H");
810 CHKERR add_parent_field(fe, UDO::OPCOL,
"G");
813 CHKERR add_parent_field(fe, UDO::OPROW,
"G");
816 CHKERR add_parent_field(fe, UDO::OPCOL,
"H");
821 auto create_subdm = [&]() {
822 auto level_ents_ptr = boost::make_shared<Range>();
828 CHKERR DMSetType(subdm,
"DMMOFEM");
834 for (
auto f : {
"H",
"G"}) {
840 if constexpr (
debug) {
844 dm_ents.subset_by_type(MBVERTEX));
845 dm_ents = subtract(dm_ents, dm_ents.subset_by_type(MBVERTEX));
853 auto subdm = create_subdm();
855 pip_mng->getDomainRhsFE().reset();
856 pip_mng->getDomainLhsFE().reset();
859 pip_mng->getDomainLhsFE()->exeTestHook = exe_test;
860 pip_mng->getDomainRhsFE()->exeTestHook = exe_test;
862 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
863 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
866 auto snes = pip_mng->createSNES(subdm);
869 auto set_section_monitor = [&](
auto solver) {
871 CHKERR SNESMonitorSet(solver,
874 (
void *)(snes_ctx_ptr.get()),
nullptr);
879 auto print_section_field = [&]() {
885 CHKERR PetscSectionGetNumFields(section, &num_fields);
886 for (
int f = 0;
f < num_fields; ++
f) {
890 <<
"Field " <<
f <<
" " << std::string(
field_name);
896 CHKERR set_section_monitor(snes);
897 CHKERR print_section_field();
899 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
903 CHKERR SNESSetFromOptions(snes);
904 CHKERR SNESSolve(snes, PETSC_NULLPTR,
D);
906 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
907 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
918 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
929 pip_mng->getOpDomainRhsPipeline().clear();
930 pip_mng->getOpDomainLhsPipeline().clear();
933 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM_X",
935 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM_X",
937 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM_Y",
939 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM_Y",
941 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX",
"U",
943 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX",
"L",
945 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"ZERO",
967 auto fe_domain_lhs = pip_mng->getDomainLhsFE();
968 auto fe_domain_rhs = pip_mng->getDomainRhsFE();
969 auto fe_bdy_lhs = pip_mng->getBoundaryLhsFE();
970 auto fe_bdy_rhs = pip_mng->getBoundaryRhsFE();
972 pip_mng->getDomainLhsFE().reset();
973 pip_mng->getDomainRhsFE().reset();
974 pip_mng->getBoundaryLhsFE().reset();
975 pip_mng->getBoundaryRhsFE().reset();
977 using UDO = ForcesAndSourcesCore::UserDataOperator;
980 auto add_parent_field_domain = [&](
auto fe,
auto op,
auto field,
auto bit) {
985 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
994 auto add_parent_field_bdy = [&](
auto fe,
auto op,
auto field,
auto bit) {
999 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1009 auto create_run_parent_op = [&](
auto parent_fe_ptr,
auto this_fe_ptr,
1011 auto parent_mask = fe_bit;
1019 auto get_parents_vel_fe_ptr = [&](
auto this_fe_ptr,
auto fe_bit) {
1020 std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
1022 parents_elems_ptr_vec.emplace_back(
1023 boost::make_shared<DomianParentEle>(
mField));
1025 parents_elems_ptr_vec[
l - 1]->getOpPtrVector().push_back(
1026 create_run_parent_op(parents_elems_ptr_vec[
l], this_fe_ptr, fe_bit));
1028 return parents_elems_ptr_vec[0];
1032 auto solve_projection = [&](
auto exe_test) {
1035 auto set_domain_rhs = [&](
auto fe) {
1037 auto &pip = fe->getOpPtrVector();
1039 auto u_ptr = boost::make_shared<MatrixDouble>();
1040 auto p_ptr = boost::make_shared<VectorDouble>();
1041 auto h_ptr = boost::make_shared<VectorDouble>();
1042 auto g_ptr = boost::make_shared<VectorDouble>();
1044 auto eval_fe_ptr = boost::make_shared<DomianParentEle>(
mField);
1046 auto &pip = eval_fe_ptr->getOpPtrVector();
1052 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1060 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW,
"U",
1063 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW,
"P",
1066 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW,
"H",
1069 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW,
"G",
1073 auto parent_eval_fe_ptr =
1075 pip.push_back(create_run_parent_op(parent_eval_fe_ptr, eval_fe_ptr,
1078 auto assemble_fe_ptr = boost::make_shared<DomianParentEle>(
mField);
1080 auto &pip = assemble_fe_ptr->getOpPtrVector();
1086 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1092 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW,
"U",
1095 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW,
"P",
1098 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW,
"H",
1101 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW,
"G",
1105 auto parent_assemble_fe_ptr =
1107 pip.push_back(create_run_parent_op(
1113 auto set_domain_lhs = [&](
auto fe) {
1116 auto &pip = fe->getOpPtrVector();
1124 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1133 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"U",
1135 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"U",
1138 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"P",
1140 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"P",
1143 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"H",
1145 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"H",
1148 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"G",
1150 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"G",
1157 auto set_bdy_rhs = [&](
auto fe) {
1159 auto &pip = fe->getOpPtrVector();
1161 auto l_ptr = boost::make_shared<VectorDouble>();
1163 auto eval_fe_ptr = boost::make_shared<BoundaryParentEle>(
mField);
1165 auto &pip = eval_fe_ptr->getOpPtrVector();
1170 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1178 CHKERR add_parent_field_bdy(eval_fe_ptr, UDO::OPROW,
"L",
1182 auto parent_eval_fe_ptr =
1184 pip.push_back(create_run_parent_op(parent_eval_fe_ptr, eval_fe_ptr,
1187 auto assemble_fe_ptr = boost::make_shared<BoundaryParentEle>(
mField);
1189 auto &pip = assemble_fe_ptr->getOpPtrVector();
1194 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1202 OpLSize(boost::shared_ptr<VectorDouble> l_ptr)
1206 if (lPtr->size() != getGaussPts().size2()) {
1207 lPtr->resize(getGaussPts().size2());
1214 boost::shared_ptr<VectorDouble> lPtr;
1217 pip.push_back(
new OpLSize(l_ptr));
1219 CHKERR add_parent_field_bdy(assemble_fe_ptr, UDO::OPROW,
"L",
1223 auto parent_assemble_fe_ptr =
1225 pip.push_back(create_run_parent_op(
1231 auto set_bdy_lhs = [&](
auto fe) {
1234 auto &pip = fe->getOpPtrVector();
1240 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1249 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"L",
1251 CHKERR add_parent_field_bdy(fe, UDO::OPCOL,
"L",
1259 auto level_ents_ptr = boost::make_shared<Range>();
1263 auto get_prj_ents = [&]() {
1268 auto common_ents = intersect(prj_mesh, *level_ents_ptr);
1269 prj_mesh = subtract(unite(*level_ents_ptr, prj_mesh), common_ents)
1275 auto prj_ents = get_prj_ents();
1279 auto rebuild = [&]() {
1283 std::vector<std::string> fields{
"U",
"P",
"H",
"G",
"L"};
1284 std::map<std::string, boost::shared_ptr<Range>> range_maps{
1286 {
"U", level_ents_ptr},
1287 {
"P", level_ents_ptr},
1288 {
"H", level_ents_ptr},
1289 {
"G", level_ents_ptr},
1290 {
"L", level_ents_ptr}
1295 simple->getProblemName(), PETSC_TRUE,
1296 &range_maps, &range_maps);
1299 CHKERR prb_mng->partitionFiniteElements(
"SUB_SOLVER",
true, 0,
1302 CHKERR prb_mng->partitionGhostDofsOnDistributedMesh(
"SUB_SOLVER");
1307 MOFEM_LOG(
"FS", Sev::verbose) <<
"Create projection problem";
1311 auto dm =
simple->getDM();
1314 CHKERR DMSetType(subdm,
"DMMOFEM");
1319 MOFEM_LOG(
"FS", Sev::inform) <<
"Nothing to project";
1324 auto create_dummy_dm = [&]() {
1327 simple->getProblemName().c_str(),
1333 auto subdm = create_subdm();
1336 pip_mng->getDomainRhsFE().reset();
1337 pip_mng->getDomainLhsFE().reset();
1338 pip_mng->getBoundaryRhsFE().reset();
1339 pip_mng->getBoundaryLhsFE().reset();
1344 pip_mng->getDomainLhsFE()->exeTestHook = exe_test;
1345 pip_mng->getDomainRhsFE()->exeTestHook = [](
FEMethod *fe_ptr) {
1348 pip_mng->getBoundaryLhsFE()->exeTestHook = exe_test;
1349 pip_mng->getBoundaryRhsFE()->exeTestHook = [](
FEMethod *fe_ptr) {
1353 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
1354 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
1355 CHKERR set_bdy_rhs(pip_mng->getCastBoundaryRhsFE());
1356 CHKERR set_bdy_lhs(pip_mng->getCastBoundaryLhsFE());
1361 auto ksp = pip_mng->createKSP(subdm);
1362 CHKERR KSPSetFromOptions(ksp);
1366 auto get_norm = [&](
auto x) {
1368 CHKERR VecNorm(x, NORM_2, &nrm);
1375 auto zero_dofs = [](boost::shared_ptr<FieldEntity> ent_ptr) {
1377 for (
auto &
v : ent_ptr->getEntFieldData()) {
1383 auto solve = [&](
auto S) {
1385 CHKERR VecZeroEntries(S);
1387 CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
1388 CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
1390 CHKERR VecGhostUpdateBegin(S, INSERT_VALUES, SCATTER_FORWARD);
1391 CHKERR VecGhostUpdateEnd(S, INSERT_VALUES, SCATTER_FORWARD);
1398 MOFEM_LOG(
"FS", Sev::inform) <<
"Solve projection";
1403 auto dummy_dm = create_dummy_dm();
1408 auto apply_restrict = [&]() {
1409 auto get_is = [](
auto v) {
1411 auto create = [&]() {
1415 CHKERR VecGetOwnershipRange(
v, &ystart, NULL);
1416 CHKERR ISCreateStride(PETSC_COMM_SELF,
n, ystart, 1, &iy);
1423 auto iy = get_is(glob_x);
1427 DMSubDomainRestrict(
simple->getDM(), s, PETSC_NULLPTR, dummy_dm),
1433 DMGetNamedGlobalVector(dummy_dm,
"TSTheta_Xdot", &Xdot),
1436 auto forward_ghost = [](
auto g) {
1438 CHKERR VecGhostUpdateBegin(
g, INSERT_VALUES, SCATTER_FORWARD);
1439 CHKERR VecGhostUpdateEnd(
g, INSERT_VALUES, SCATTER_FORWARD);
1446 if constexpr (
debug) {
1448 <<
"Reverse restrict: X0 " << get_norm(X0) <<
" Xdot "
1452 return std::vector<Vec>{X0, Xdot};
1455 auto ts_solver_vecs = apply_restrict();
1457 if (ts_solver_vecs.size()) {
1459 for (
auto v : ts_solver_vecs) {
1460 MOFEM_LOG(
"FS", Sev::inform) <<
"Solve projection vector";
1466 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
1467 MOFEM_LOG(
"WORLD", Sev::verbose) <<
"Zero field " << f;
1468 CHKERR field_blas->fieldLambdaOnEntities(zero_dofs, f);
1476 MOFEM_LOG(
"FS", Sev::inform) <<
"Norm V " << get_norm(
v);
1479 CHKERR DMRestoreNamedGlobalVector(dummy_dm,
"TSTheta_X0",
1480 &ts_solver_vecs[0]);
1481 CHKERR DMRestoreNamedGlobalVector(dummy_dm,
"TSTheta_Xdot",
1482 &ts_solver_vecs[1]);
1485 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
1486 MOFEM_LOG(
"WORLD", Sev::verbose) <<
"Zero field " << f;
1487 CHKERR field_blas->fieldLambdaOnEntities(zero_dofs, f);
1496 auto post_proc = [&](
auto exe_test) {
1498 auto post_proc_fe = boost::make_shared<PostProcEleDomain>(
mField);
1499 auto &pip = post_proc_fe->getOpPtrVector();
1500 post_proc_fe->exeTestHook = exe_test;
1508 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1511 fe_parent->getOpPtrVector(), {H1});
1519 auto u_ptr = boost::make_shared<MatrixDouble>();
1520 auto p_ptr = boost::make_shared<VectorDouble>();
1521 auto h_ptr = boost::make_shared<VectorDouble>();
1522 auto g_ptr = boost::make_shared<VectorDouble>();
1524 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"U",
1527 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"P",
1530 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"H",
1533 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"G",
1537 post_proc_fe->getOpPtrVector().push_back(
1539 new OpPPMap(post_proc_fe->getPostProcMesh(),
1540 post_proc_fe->getMapGaussPts(),
1542 {{
"P", p_ptr}, {
"H", h_ptr}, {
"G", g_ptr}},
1554 auto dm =
simple->getDM();
1556 CHKERR post_proc_fe->writeFile(
"out_projection.h5m");
1565 if constexpr (
debug) {
1571 fe_domain_lhs.swap(pip_mng->getDomainLhsFE());
1572 fe_domain_rhs.swap(pip_mng->getDomainRhsFE());
1573 fe_bdy_lhs.swap(pip_mng->getBoundaryLhsFE());
1574 fe_bdy_rhs.swap(pip_mng->getBoundaryRhsFE());
1585 using UDO = ForcesAndSourcesCore::UserDataOperator;
1587 auto add_parent_field_domain = [&](
auto fe,
auto op,
auto field) {
1592 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1600 auto add_parent_field_bdy = [&](
auto fe,
auto op,
auto field) {
1605 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1618 auto dot_u_ptr = boost::make_shared<MatrixDouble>();
1619 auto u_ptr = boost::make_shared<MatrixDouble>();
1620 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
1621 auto dot_h_ptr = boost::make_shared<VectorDouble>();
1622 auto h_ptr = boost::make_shared<VectorDouble>();
1623 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
1624 auto g_ptr = boost::make_shared<VectorDouble>();
1625 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
1626 auto lambda_ptr = boost::make_shared<VectorDouble>();
1627 auto p_ptr = boost::make_shared<VectorDouble>();
1628 auto div_u_ptr = boost::make_shared<VectorDouble>();
1632 auto set_domain_general = [&](
auto fe) {
1634 auto &pip = fe->getOpPtrVector();
1642 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1645 fe_parent->getOpPtrVector(), {H1});
1651 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"U");
1671 "Coordinate system not implemented");
1674 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"H");
1680 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"G");
1685 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"P");
1690 auto set_domain_rhs = [&](
auto fe) {
1692 auto &pip = fe->getOpPtrVector();
1694 CHKERR set_domain_general(fe);
1696 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"U");
1697 pip.push_back(
new OpRhsU(
"U", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr,
1698 grad_h_ptr, g_ptr, p_ptr));
1700 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"H");
1701 pip.push_back(
new OpRhsH<false>(
"H", u_ptr, dot_h_ptr, h_ptr, grad_h_ptr,
1704 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"G");
1705 pip.push_back(
new OpRhsG<false>(
"G", h_ptr, grad_h_ptr, g_ptr));
1707 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"P");
1709 "P", div_u_ptr, [](
const double r,
const double,
const double) {
1713 "P", p_ptr, [](
const double r,
const double,
const double) {
1719 auto set_domain_lhs = [&](
auto fe) {
1721 auto &pip = fe->getOpPtrVector();
1723 CHKERR set_domain_general(fe);
1725 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"U");
1727 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"U");
1728 pip.push_back(
new OpLhsU_dU(
"U", u_ptr, grad_u_ptr, h_ptr));
1729 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"H");
1731 new OpLhsU_dH(
"U",
"H", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr, g_ptr));
1733 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"G");
1734 pip.push_back(
new OpLhsU_dG(
"U",
"G", grad_h_ptr));
1737 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"H");
1739 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"U");
1740 pip.push_back(
new OpLhsH_dU(
"H",
"U", grad_h_ptr));
1741 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"H");
1743 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"G");
1747 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"G");
1749 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"H");
1751 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"G");
1755 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"P");
1757 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"U");
1763 [](
const double r,
const double,
const double) {
1771 [](
const double r,
const double,
const double) {
1778 "Coordinate system not implemented");
1781 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"P");
1782 pip.push_back(
new OpDomainMassP(
"P",
"P", [](
double r,
double,
double) {
1790 auto get_block_name = [](
auto name) {
1791 return boost::format(
"%s(.*)") %
"WETTING_ANGLE";
1794 auto get_blocks = [&](
auto &&name) {
1796 std::regex(name.str()));
1799 auto set_boundary_rhs = [&](
auto fe) {
1801 auto &pip = fe->getOpPtrVector();
1807 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1814 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"U");
1817 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"L");
1822 pip,
mField,
"L", {},
"INFLUX",
1825 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"U");
1828 auto wetting_block = get_blocks(get_block_name(
"WETTING_ANGLE"));
1829 if (wetting_block.size()) {
1836 op_bdy_side->getOpPtrVector(), {H1});
1839 op_bdy_side->getSideFEPtr(),
"", UDO::OPSPACE,
1843 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1846 fe_parent->getOpPtrVector(), {H1});
1852 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPROW,
1854 op_bdy_side->getOpPtrVector().push_back(
1857 pip.push_back(op_bdy_side);
1860 for (
auto &b : wetting_block) {
1862 std::vector<double> attr_vec;
1863 CHKERR b->getMeshsetIdEntitiesByDimension(
1865 b->getAttributes(attr_vec);
1866 if (attr_vec.size() != 1)
1868 "Should be one attribute");
1869 MOFEM_LOG(
"FS", Sev::inform) <<
"Wetting angle: " << attr_vec.front();
1871 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"G");
1873 "G", grad_h_ptr, boost::make_shared<Range>(force_edges),
1881 auto set_boundary_lhs = [&](
auto fe) {
1883 auto &pip = fe->getOpPtrVector();
1889 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1896 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"L");
1897 CHKERR add_parent_field_bdy(fe, UDO::OPCOL,
"U");
1900 auto wetting_block = get_blocks(get_block_name(
"WETTING_ANGLE"));
1901 if (wetting_block.size()) {
1902 auto col_ind_ptr = boost::make_shared<std::vector<VectorInt>>();
1903 auto col_diff_base_ptr = boost::make_shared<std::vector<MatrixDouble>>();
1911 op_bdy_side->getOpPtrVector(), {H1});
1914 op_bdy_side->getSideFEPtr(),
"", UDO::OPSPACE,
1918 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1921 fe_parent->getOpPtrVector(), {H1});
1927 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPROW,
1929 op_bdy_side->getOpPtrVector().push_back(
1931 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPCOL,
1933 op_bdy_side->getOpPtrVector().push_back(
1937 pip.push_back(op_bdy_side);
1940 for (
auto &b : wetting_block) {
1942 std::vector<double> attr_vec;
1943 CHKERR b->getMeshsetIdEntitiesByDimension(
1945 b->getAttributes(attr_vec);
1946 if (attr_vec.size() != 1)
1948 "Should be one attribute");
1950 <<
"wetting angle edges size " << force_edges.size();
1952 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"G");
1954 "G", grad_h_ptr, col_ind_ptr, col_diff_base_ptr,
1955 boost::make_shared<Range>(force_edges), attr_vec.front()));
1964 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
1965 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
1966 CHKERR set_boundary_rhs(pip_mng->getCastBoundaryRhsFE());
1967 CHKERR set_boundary_lhs(pip_mng->getCastBoundaryLhsFE());
1991 boost::shared_ptr<PostProcEleDomainCont> post_proc,
1992 boost::shared_ptr<PostProcEleBdyCont> post_proc_edge,
1993 std::pair<boost::shared_ptr<BoundaryEle>, boost::shared_ptr<VectorDouble>>
2000 MOFEM_LOG(
"FS", Sev::verbose) <<
"Monitor";
2003 MOFEM_LOG(
"FS", Sev::verbose) <<
"Mesh pre proc";
2006 auto post_proc_begin =
2007 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
2009 auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
2017 CHKERR post_proc_end->writeFile(
2018 "out_step_" + boost::lexical_cast<std::string>(
ts_step) +
".h5m");
2019 MOFEM_LOG(
"FS", Sev::verbose) <<
"Mesh pre proc done";
2025 MPI_Allreduce(MPI_IN_PLACE, &(*
liftVec)[0],
SPACE_DIM, MPI_DOUBLE, MPI_SUM,
2029 <<
" lift vec x: " << (*liftVec)[0] <<
" y: " << (*liftVec)[1];
2047 using UDO = ForcesAndSourcesCore::UserDataOperator;
2055 auto setup_subdm = [&](
auto dm) {
2059 auto dm =
simple->getDM();
2060 auto level_ents_ptr = boost::make_shared<Range>();
2061 CHKERR bit_mng->getEntitiesByRefLevel(
2064 CHKERR DMSetType(subdm,
"DMMOFEM");
2070 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
2083 auto get_fe_post_proc = [&](
auto post_proc_mesh) {
2084 auto add_parent_field_domain = [&](
auto fe,
auto op,
auto field) {
2089 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2098 boost::make_shared<PostProcEleDomainCont>(
mField, post_proc_mesh);
2099 post_proc_fe->exeTestHook = [](
FEMethod *fe_ptr) {
2104 auto u_ptr = boost::make_shared<MatrixDouble>();
2105 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
2106 auto h_ptr = boost::make_shared<VectorDouble>();
2107 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
2108 auto p_ptr = boost::make_shared<VectorDouble>();
2109 auto g_ptr = boost::make_shared<VectorDouble>();
2110 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
2113 post_proc_fe->getOpPtrVector(), {H1});
2119 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2122 fe_parent->getOpPtrVector(), {H1});
2128 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"U");
2129 post_proc_fe->getOpPtrVector().push_back(
2131 post_proc_fe->getOpPtrVector().push_back(
2135 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"H");
2136 post_proc_fe->getOpPtrVector().push_back(
2138 post_proc_fe->getOpPtrVector().push_back(
2141 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"P");
2142 post_proc_fe->getOpPtrVector().push_back(
2145 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"G");
2146 post_proc_fe->getOpPtrVector().push_back(
2148 post_proc_fe->getOpPtrVector().push_back(
2153 post_proc_fe->getOpPtrVector().push_back(
2156 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
2158 {{
"H", h_ptr}, {
"P", p_ptr}, {
"G", g_ptr}},
2160 {{
"U", u_ptr}, {
"H_GRAD", grad_h_ptr}, {
"G_GRAD", grad_g_ptr}},
2162 {{
"GRAD_U", grad_u_ptr}},
2170 return post_proc_fe;
2173 auto get_bdy_post_proc_fe = [&](
auto post_proc_mesh) {
2174 auto add_parent_field_bdy = [&](
auto fe,
auto op,
auto field) {
2175 return setParentDofs(
2179 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2188 boost::make_shared<PostProcEleBdyCont>(mField, post_proc_mesh);
2189 post_proc_fe->exeTestHook = [](
FEMethod *fe_ptr) {
2198 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2207 OpGetNormal(boost::shared_ptr<VectorDouble> l_ptr,
2208 boost::shared_ptr<MatrixDouble> n_ptr)
2216 auto t_n_fe = getFTensor1NormalsAtGaussPts();
2217 ptrNormal->resize(
SPACE_DIM, getGaussPts().size2(),
false);
2219 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
2220 t_n(
i) = t_n_fe(
i) * t_l / std::sqrt(t_n_fe(
i) * t_n_fe(
i));
2229 boost::shared_ptr<VectorDouble> ptrL;
2230 boost::shared_ptr<MatrixDouble> ptrNormal;
2233 auto u_ptr = boost::make_shared<MatrixDouble>();
2234 auto p_ptr = boost::make_shared<VectorDouble>();
2235 auto lambda_ptr = boost::make_shared<VectorDouble>();
2236 auto normal_l_ptr = boost::make_shared<MatrixDouble>();
2238 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW,
"U");
2239 post_proc_fe->getOpPtrVector().push_back(
2242 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW,
"L");
2243 post_proc_fe->getOpPtrVector().push_back(
2245 post_proc_fe->getOpPtrVector().push_back(
2246 new OpGetNormal(lambda_ptr, normal_l_ptr));
2248 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW,
"P");
2249 post_proc_fe->getOpPtrVector().push_back(
2253 op_ptr->doWorkRhsHook = [&](
DataOperator *base_op_ptr,
int side,
2262 post_proc_fe->getOpPtrVector().push_back(
2264 new OpPPMap(post_proc_fe->getPostProcMesh(),
2265 post_proc_fe->getMapGaussPts(),
2279 return post_proc_fe;
2282 auto get_lift_fe = [&]() {
2283 auto add_parent_field_bdy = [&](
auto fe,
auto op,
auto field) {
2284 return setParentDofs(
2288 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2296 auto fe = boost::make_shared<BoundaryEle>(mField);
2297 fe->exeTestHook = [](
FEMethod *fe_ptr) {
2302 auto lift_ptr = boost::make_shared<VectorDouble>();
2303 auto p_ptr = boost::make_shared<VectorDouble>();
2304 auto ents_ptr = boost::make_shared<Range>();
2310 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2317 std::vector<const CubitMeshSets *> vec_ptr;
2319 std::regex(
"LIFT"), vec_ptr);
2320 for (
auto m_ptr : vec_ptr) {
2321 auto meshset = m_ptr->getMeshset();
2325 ents_ptr->merge(ents);
2328 MOFEM_LOG(
"FS", Sev::noisy) <<
"Lift ents " << (*ents_ptr);
2330 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"P");
2331 fe->getOpPtrVector().push_back(
2333 fe->getOpPtrVector().push_back(
2336 return std::make_pair(fe, lift_ptr);
2339 auto set_post_proc_monitor = [&](
auto ts) {
2343 boost::shared_ptr<FEMethod> null_fe;
2344 auto post_proc_mesh = boost::make_shared<moab::Core>();
2345 auto monitor_ptr = boost::make_shared<Monitor>(
2347 get_fe_post_proc(post_proc_mesh), get_bdy_post_proc_fe(post_proc_mesh),
2350 null_fe, monitor_ptr);
2354 auto dm =
simple->getDM();
2358 auto ts_pre_post_proc = boost::make_shared<TSPrePostProc>();
2363 ptr->fsRawPtr =
this;
2364 ptr->solverSubDM = create_solver_dm(
simple->getDM());
2368 CHKERR VecAssemblyBegin(ptr->globSol);
2369 CHKERR VecAssemblyEnd(ptr->globSol);
2371 auto sub_ts = pip_mng->createTSIM(ptr->solverSubDM);
2373 CHKERR set_post_proc_monitor(sub_ts);
2376 CHKERR TSSetFromOptions(ts);
2380 auto print_fields_in_section = [&]() {
2383 simple->getProblemName());
2384 PetscInt num_fields;
2385 CHKERR PetscSectionGetNumFields(section, &num_fields);
2386 for (
int f = 0;
f < num_fields; ++
f) {
2390 <<
"Field " <<
f <<
" " << std::string(
field_name);
2395 CHKERR print_fields_in_section();
2397 CHKERR TSSolve(ts, ptr->globSol);
2405#ifdef PYTHON_INIT_SURFACE
2410 const char param_file[] =
"param_file.petsc";
2414 auto core_log = logging::core::get();
2422 DMType dm_name =
"DMMOFEM";
2427 moab::Core mb_instance;
2428 moab::Interface &moab = mb_instance;
2445#ifdef PYTHON_INIT_SURFACE
2446 if (Py_FinalizeEx() < 0) {
2462 "can not get vertices on bit 0");
2467 Range plus_range, minus_range;
2468 std::vector<EntityHandle> plus, minus;
2471 for (
auto p = vertices.pair_begin(); p != vertices.pair_end(); ++p) {
2473 const auto f = p->first;
2474 const auto s = p->second;
2479 auto it = dofs_mi.lower_bound(lo_uid);
2480 const auto hi_it = dofs_mi.upper_bound(hi_uid);
2484 plus.reserve(std::distance(it, hi_it));
2485 minus.reserve(std::distance(it, hi_it));
2487 for (; it != hi_it; ++it) {
2488 const auto v = (*it)->getFieldData();
2490 plus.push_back((*it)->getEnt());
2492 minus.push_back((*it)->getEnt());
2495 plus_range.insert_list(plus.begin(), plus.end());
2496 minus_range.insert_list(minus.begin(), minus.end());
2501 <<
"Plus range " << plus_range << endl;
2503 <<
"Minus range " << minus_range << endl;
2506 auto get_elems = [&](
auto &ents,
auto bit,
auto mask) {
2509 moab::Interface::UNION),
2510 "can not get adjacencies");
2512 "can not filter elements with bit 0");
2516 CHKERR comm_mng->synchroniseEntities(plus_range);
2517 CHKERR comm_mng->synchroniseEntities(minus_range);
2520 ele_plus[0] = get_elems(plus_range,
bit(0),
BitRefLevel().set());
2521 ele_minus[0] = get_elems(minus_range,
bit(0),
BitRefLevel().set());
2522 auto common = intersect(ele_plus[0], ele_minus[0]);
2523 ele_plus[0] = subtract(ele_plus[0], common);
2524 ele_minus[0] = subtract(ele_minus[0], common);
2526 auto get_children = [&](
auto &p,
auto &
c) {
2528 CHKERR bit_mng->updateRangeByChildren(p,
c);
2540 auto get_level = [&](
auto &p,
auto &
m,
auto z,
auto bit,
auto mask) {
2543 bit_mng->getEntitiesByDimAndRefLevel(
bit, mask,
SPACE_DIM,
l),
2544 "can not get vertices on bit");
2547 for (
auto f = 0; f != z; ++f) {
2551 l = get_elems(conn,
bit, mask);
2556 std::vector<Range> vec_levels(
nb_levels);
2558 vec_levels[
l] = get_level(ele_plus[
l], ele_minus[
l], 2 * overlap,
bit(
l),
2562 if constexpr (
debug) {
2565 std::string name = (boost::format(
"out_r%d.h5m") %
l).str();
2582 start_mask[s] =
true;
2588 Range prev_level_skin;
2592 CHKERR bit_mng->lambdaBitRefLevel(
2599 auto set_levels = [&](
auto &&
2609 auto get_adj = [&](
auto ents) {
2615 ents.subset_by_dimension(
SPACE_DIM), d,
false, conn,
2616 moab::Interface::UNION),
2630 CHKERR bit_mng->updateRangeByParent(vec_levels[
l], parents);
2633 level_prev = subtract(level_prev, parents);
2635 level_prev.merge(vec_levels[
l]);
2643 CHKERR bit_mng->getEntitiesByDimAndRefLevel(
2645 level = get_adj(level);
2654 auto set_skins = [&]() {
2658 ParallelComm *pcomm =
2667 "can't get bit level");
2674 auto get_level_skin = [&]() {
2680 CHKERR pcomm->filter_pstatus(skin_level_mesh,
2681 PSTATUS_SHARED | PSTATUS_MULTISHARED,
2682 PSTATUS_NOT, -1,
nullptr);
2685 skin_level = subtract(skin_level,
2689 Range skin_level_verts;
2692 skin_level.merge(skin_level_verts);
2695 bit_prev.set(
l - 1);
2699 skin.merge(subtract(skin_level, level_prev));
2705 auto resolve_shared = [&](
auto &&skin) {
2706 Range tmp_skin = skin;
2708 map<int, Range> map_procs;
2710 tmp_skin, &map_procs);
2712 Range from_other_procs;
2713 for (
auto &
m : map_procs) {
2715 from_other_procs.merge(
m.second);
2719 auto common = intersect(
2720 skin, from_other_procs);
2721 skin.merge(from_other_procs);
2725 if (!common.empty()) {
2727 skin = subtract(skin, common);
2734 auto get_parent_level_skin = [&](
auto skin) {
2736 CHKERR bit_mng->updateRangeByParent(
2737 skin.subset_by_dimension(
SPACE_DIM - 1), skin_parents);
2738 Range skin_parent_verts;
2741 skin_parents.merge(skin_parent_verts);
2744 return skin_parents;
2747 auto child_skin = resolve_shared(get_level_skin());
2748 auto parent_skin = get_parent_level_skin(child_skin);
2750 child_skin = subtract(child_skin, parent_skin);
2758 auto set_current = [&]() {
2767 last_level = subtract(last_level, skin_child);
2779 if constexpr (
debug) {
2782 std::string name = (boost::format(
"out_level%d.h5m") %
l).str();
2785 "PARALLEL=WRITE_PART");
2789 "MOAB",
"PARALLEL=WRITE_PART");
2792 "MOAB",
"PARALLEL=WRITE_PART");
2796 "MOAB",
"PARALLEL=WRITE_PART");
2799 "MOAB",
"PARALLEL=WRITE_PART");
2807 boost::shared_ptr<FEMethod> fe_top, std::string
field_name,
2808 ForcesAndSourcesCore::UserDataOperator::OpType op,
2810 boost::function<boost::shared_ptr<ForcesAndSourcesCore>()> get_elem,
2817 boost::function<void(boost::shared_ptr<ForcesAndSourcesCore>,
int)>
2819 [&](boost::shared_ptr<ForcesAndSourcesCore> parent_fe_pt,
int level) {
2824 auto fe_ptr_current = get_elem();
2828 boost::dynamic_pointer_cast<ForcesAndSourcesCore>(
2833 if (op == ForcesAndSourcesCore::UserDataOperator::OPSPACE) {
2836 parent_fe_pt->getOpPtrVector().push_back(
2840 H1, op, fe_ptr_current,
2851 parent_fe_pt->getOpPtrVector().push_back(
2866 add_parent_level(boost::dynamic_pointer_cast<ForcesAndSourcesCore>(fe_top),
2881 auto cut_off_dofs = [&]() {
2884 auto &m_field = ptr->fsRawPtr->mField;
2886 Range current_verts;
2891 auto cut_off_verts = [&](boost::shared_ptr<FieldEntity> ent_ptr) {
2893 for (
auto &
h : ent_ptr->getEntFieldData()) {
2899 auto field_blas = m_field.getInterface<
FieldBlas>();
2909 MOFEM_LOG(
"FS", Sev::inform) <<
"Run step pre proc";
2911 auto &m_field = ptr->fsRawPtr->mField;
2915 auto get_norm = [&](
auto x) {
2917 CHKERR VecNorm(x, NORM_2, &nrm);
2922 auto refine_problem = [&]() {
2924 MOFEM_LOG(
"FS", Sev::inform) <<
"Refine problem";
2926 CHKERR ptr->fsRawPtr->projectData();
2932 auto set_jacobian_operators = [&]() {
2935 CHKERR KSPReset(ptr->subKSP);
2940 auto set_solution = [&]() {
2942 MOFEM_LOG(
"FS", Sev::inform) <<
"Set solution";
2944 PetscObjectState state;
2955 INSERT_VALUES, SCATTER_FORWARD);
2958 <<
"Set solution, vector norm " << get_norm(ptr->globSol);
2963 PetscObjectTypeCompare((PetscObject)ts, TSTHETA, &is_theta);
2967 CHKERR set_jacobian_operators();
2972 "Sorry, only TSTheta handling is implemented");
2977 PetscBarrier((PetscObject)ts);
2989 auto &m_field = ptr->fsRawPtr->mField;
3001 auto sub_u = ptr->getSubVector();
3004 auto scatter = ptr->getScatter(sub_u, u,
R);
3006 auto apply_scatter_and_update = [&](
auto x,
auto sub_x) {
3008 CHKERR VecScatterBegin(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3009 CHKERR VecScatterEnd(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3010 CHKERR VecGhostUpdateBegin(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3011 CHKERR VecGhostUpdateEnd(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3015 CHKERR apply_scatter_and_update(u, sub_u);
3016 CHKERR apply_scatter_and_update(u_t, sub_u_t);
3019 CHKERR VecScatterBegin(scatter, sub_f, f, INSERT_VALUES, SCATTER_FORWARD);
3020 CHKERR VecScatterEnd(scatter, sub_f, f, INSERT_VALUES, SCATTER_FORWARD);
3026 PetscReal
a, Mat A, Mat
B,
3030 auto sub_u = ptr->getSubVector();
3032 auto scatter = ptr->getScatter(sub_u, u,
R);
3034 auto apply_scatter_and_update = [&](
auto x,
auto sub_x) {
3036 CHKERR VecScatterBegin(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3037 CHKERR VecScatterEnd(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3038 CHKERR VecGhostUpdateBegin(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3039 CHKERR VecGhostUpdateEnd(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3043 CHKERR apply_scatter_and_update(u, sub_u);
3044 CHKERR apply_scatter_and_update(u_t, sub_u_t);
3047 ptr->tsCtxPtr.get());
3056 auto get_norm = [&](
auto x) {
3058 CHKERR VecNorm(x, NORM_2, &nrm);
3062 auto sub_u = ptr->getSubVector();
3063 auto scatter = ptr->getScatter(sub_u, u,
R);
3064 CHKERR VecScatterBegin(scatter, u, sub_u, INSERT_VALUES, SCATTER_REVERSE);
3065 CHKERR VecScatterEnd(scatter, u, sub_u, INSERT_VALUES, SCATTER_REVERSE);
3066 CHKERR VecGhostUpdateBegin(sub_u, INSERT_VALUES, SCATTER_FORWARD);
3067 CHKERR VecGhostUpdateEnd(sub_u, INSERT_VALUES, SCATTER_FORWARD);
3070 <<
"u norm " << get_norm(u) <<
" u sub nom " << get_norm(sub_u);
3080 MOFEM_LOG(
"FS", Sev::verbose) <<
"SetUP sub PC";
3081 ptr->subKSP =
createKSP(ptr->fsRawPtr->mField.get_comm());
3082 CHKERR KSPSetFromOptions(ptr->subKSP);
3090 auto sub_x = ptr->getSubVector();
3092 auto scatter = ptr->getScatter(sub_x, pc_x,
R);
3093 CHKERR VecScatterBegin(scatter, pc_f, sub_f, INSERT_VALUES,
3095 CHKERR VecScatterEnd(scatter, pc_f, sub_f, INSERT_VALUES, SCATTER_REVERSE);
3096 CHKERR KSPSetOperators(ptr->subKSP, ptr->subB, ptr->subB);
3097 MOFEM_LOG(
"FS", Sev::verbose) <<
"PCShell solve";
3098 CHKERR KSPSolve(ptr->subKSP, sub_f, sub_x);
3099 MOFEM_LOG(
"FS", Sev::verbose) <<
"PCShell solve <- done";
3100 CHKERR VecScatterBegin(scatter, sub_x, pc_x, INSERT_VALUES,
3102 CHKERR VecScatterEnd(scatter, sub_x, pc_x, INSERT_VALUES, SCATTER_FORWARD);
3109 auto prb_ptr = ptr->fsRawPtr->mField.get_problem(
"SUB_SOLVER");
3110 if (
auto sub_data = prb_ptr->getSubData()) {
3111 auto is = sub_data->getSmartColIs();
3138 auto dm =
simple->getDM();
3146 CHKERR TSGetSNES(ts, &snes);
3149 auto set_section_monitor = [&](
auto snes) {
3151 CHKERR SNESMonitorSet(snes,
3154 (
void *)(snes_ctx_ptr.get()),
nullptr);
3158 CHKERR set_section_monitor(snes);
3160 auto ksp =
createKSP(m_field.get_comm());
3161 CHKERR KSPSetType(ksp, KSPPREONLY);
3163 CHKERR PCSetType(sub_pc, PCSHELL);
3166 CHKERR KSPSetPC(ksp, sub_pc);
3167 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 >::BoundaryParentEle BoundaryParentEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
[Define dimension]
BoundaryEle::UserDataOperator BoundaryEleOp
Kronecker Delta class symmetric.
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
#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
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, U_FIELD_DIM > OpDomainMassU
constexpr int U_FIELD_DIM
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
static boost::weak_ptr< TSPrePostProc > tsPrePostProc
FTensor::Index< 'j', SPACE_DIM > j
FormsIntegrators< BoundaryEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, 1 > OpBoundaryMassL
double lambda
surface tension
FTensor::Index< 'k', SPACE_DIM > k
OpDomainMassH OpDomainMassG
FTensor::Index< 'i', SPACE_DIM > i
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesVector< BASE_DIM, SPACE_DIM, 1 > OpDomainAssembleVector
auto init_h
Initialisation function.
FTensor::Index< 'l', SPACE_DIM > l
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpDomainAssembleScalar
OpDomainMassH OpDomainMassP
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMixScalarTimesDiv< SPACE_DIM, COORD_TYPE > OpMixScalarTimesDiv
auto get_dofs_ents_all
get entities of dofs in the problem - used for debugging
constexpr IntegrationType I
auto get_skin_projection_bit
auto get_dofs_ents_by_field_name
get entities of dofs in the problem - used for debugging
ElementsAndOps< SPACE_DIM >::DomianParentEle DomianParentEle
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, 1 > OpDomainMassH
auto wetting_angle_sub_stepping
int order
approximation order
SideEle::UserDataOperator SideOp
auto get_current_bit
dofs bit used to do calculations
FormsIntegrators< BoundaryEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpBoundaryAssembleScalar
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 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
PipelineManager::ElementsAndOpsByDim< FE_DIM >::DomianParentEle DomianParentEle
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)
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
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.
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.
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
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
FTensor::Index< 'm', 3 > m
MoFEM::Interface & mField
Range findParentsToRefine(Range ents, BitRefLevel level, BitRefLevel mask)
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode runProblem()
[Run programme]
FreeSurface(MoFEM::Interface &m_field)
MoFEMErrorCode projectData()
[Boundary condition]
MoFEMErrorCode boundaryCondition()
[Set up problem]
MoFEMErrorCode readMesh()
[Run programme]
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 parents levels.
std::vector< Range > findEntitiesCrossedByPhaseInterface(size_t overlap)
Find entities on refinement levels.
MoFEMErrorCode assembleSystem()
[Data projection]
MoFEM::Interface & mField
MoFEMErrorCode refineMesh(size_t overlap)
MoFEMErrorCode makeRefProblem()
MoFEMErrorCode solveSystem()
[Solve]
Add operators pushing bases from local to physical configuration.
boost::weak_ptr< CacheTuple > getCacheWeakPtr() const
Get the cache weak ptr object.
Simple interface for fast problem set-up.
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
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 field values (template specialization) for tensor field rank 1, i.e. vector field.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Get value at integration points for scalar field.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Approximate field values for given petsc vector.
Get values at integration pts for tensor field rank 1, i.e. vector field.
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...
PipelineManager interface.
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
bool & getParentAdjacencies()
Get the addParentAdjacencies.
intrusive_ptr for managing petsc objects
PetscInt ts_step
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()
function is run at the end of loop
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()
static MoFEMErrorCode tsMonitor(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Wrapper for TS monitor.
static MoFEMErrorCode pcApply(PC pc, Vec pc_f, Vec pc_x)
SmartPetscObj< VecScatter > getScatter(Vec x, Vec y, enum FR fr)
SmartPetscObj< DM > solverSubDM
static MoFEMErrorCode tsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Wrapper for SNES Lhs.
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)
static MoFEMErrorCode pcSetup(PC pc)
static MoFEMErrorCode tsPostProc(TS ts)
Post process time step.
static MoFEMErrorCode tsPreStage(TS ts)
MoFEMErrorCode tsSetUp(TS ts)
Used to setup TS solver.
SmartPetscObj< KSP > subKSP
constexpr CoordinateTypes COORD_TYPE