13#include <petsc/private/petscimpl.h>
17static char help[] =
"...\n\n";
19#ifdef PYTHON_INIT_SURFACE
20#include <boost/python.hpp>
21#include <boost/python/def.hpp>
22namespace bp = boost::python;
25 SurfacePython() =
default;
26 virtual ~SurfacePython() =
default;
33 auto main_module = bp::import(
"__main__");
34 mainNamespace = main_module.attr(
"__dict__");
35 bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
37 surfaceFun = mainNamespace[
"surface"];
39 }
catch (bp::error_already_set
const &) {
49 double x,
double y,
double z,
double eta,
double &s
56 s = bp::extract<double>(surfaceFun(x, y, z,
eta));
58 }
catch (bp::error_already_set
const &) {
67 bp::object mainNamespace;
68 bp::object surfaceFun;
71static boost::weak_ptr<SurfacePython> surfacePythonWeakPtr;
83 IntegrationType::GAUSS;
95using SideOp = SideEle::UserDataOperator;
123template <CoordinateTypes COORD_TYPE>
176double tol = std::numeric_limits<float>::epsilon();
195 constexpr int sub_stepping = 16;
196 return std::min(1.,
static_cast<double>(ts_step) / sub_stepping);
201auto my_max = [](
const double x) {
return (x - 1 + std::abs(x + 1)) / 2; };
202auto my_min = [](
const double x) {
return (x + 1 - std::abs(x - 1)) / 2; };
205 if (
h >= -1 &&
h < 1)
212 return diff *
h + ave;
221auto get_f = [](
const double h) {
return 4 *
W *
h * (
h *
h - 1); };
228 return md * (1 -
h *
h);
234 const double h2 =
h *
h;
235 const double h3 = h2 *
h;
237 return md * (2 * h3 - 3 * h2 + 1);
239 return md * (-2 * h3 - 3 * h2 + 1);
244 return md * (6 *
h * (
h - 1));
246 return md * (-6 *
h * (
h + 1));
262 constexpr double R = 0.0125;
263 constexpr double A =
R * 0.2;
264 const double theta = atan2(r, y);
265 const double w =
R +
A * cos(
n * theta);
266 const double d = std::sqrt(r * r + y * y);
267 return tanh((w - d) / (
eta * std::sqrt(2)));
271 constexpr double y0 = 0.5;
272 constexpr double R = 0.4;
274 const double d = std::sqrt(r * r + y * y);
275 return tanh((
R - d) / (
eta * std::sqrt(2)));
279 constexpr double water_height = 0.;
280 return tanh((water_height - y) / (
eta * std::sqrt(2)));
285 return -tanh((-0.039 - x) / (
eta * std::sqrt(2)));
295auto init_h = [](
double r,
double y,
double theta) {
296#ifdef PYTHON_INIT_SURFACE
298 if (
auto ptr = surfacePythonWeakPtr.lock()) {
300 "error eval python");
315 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel();
320 MPI_Allreduce(&l_size, &g_size, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
324auto save_range = [](moab::Interface &moab,
const std::string name,
329 CHKERR moab.add_entities(*out_meshset, r);
330 CHKERR moab.write_file(name.c_str(),
"MOAB",
"PARALLEL=WRITE_PART",
331 out_meshset->get_ptr(), 1);
342 std::vector<EntityHandle> ents_vec;
348 auto dofs = prb_ptr->numeredRowDofsPtr;
351 ents_vec.reserve(std::distance(lo_it, hi_it));
353 for (; lo_it != hi_it; ++lo_it) {
354 ents_vec.push_back((*lo_it)->getEnt());
357 std::sort(ents_vec.begin(), ents_vec.end());
358 auto it = std::unique(ents_vec.begin(), ents_vec.end());
360 r.insert_list(ents_vec.begin(), it);
370 std::vector<EntityHandle> ents_vec;
375 auto dofs = prb_ptr->numeredRowDofsPtr;
376 ents_vec.reserve(dofs->size());
378 for (
auto d : *dofs) {
379 ents_vec.push_back(d->getEnt());
382 std::sort(ents_vec.begin(), ents_vec.end());
383 auto it = std::unique(ents_vec.begin(), ents_vec.end());
385 r.insert_list(ents_vec.begin(), it);
450 PetscReal a, Mat
A, Mat
B,
461 boost::shared_ptr<SnesCtx>
463 boost::shared_ptr<TsCtx>
512 boost::shared_ptr<FEMethod> fe_top, std::string
field_name,
513 ForcesAndSourcesCore::UserDataOperator::OpType op,
515 boost::function<boost::shared_ptr<ForcesAndSourcesCore>()> get_elem,
536 MOFEM_LOG(
"FS", Sev::inform) <<
"Read mesh for problem";
539 simple->getParentAdjacencies() =
true;
559 const char *coord_type_names[] = {
"cartesian",
"polar",
"cylindrical",
564 MOFEM_LOG(
"FS", Sev::inform) <<
"Approximation order = " <<
order;
566 <<
"Number of refinement levels nb_levels = " <<
nb_levels;
574 "-rho_m", &
rho_m, PETSC_NULL);
584 "Height of the well in energy functional",
"-W",
598 MOFEM_LOG(
"FS", Sev::inform) <<
"Acceleration a0 = " <<
a0;
599 MOFEM_LOG(
"FS", Sev::inform) <<
"\"Minus\" phase density rho_m = " <<
rho_m;
600 MOFEM_LOG(
"FS", Sev::inform) <<
"\"Minus\" phase viscosity mu_m = " <<
mu_m;
601 MOFEM_LOG(
"FS", Sev::inform) <<
"\"Plus\" phase density rho_p = " <<
rho_p;
602 MOFEM_LOG(
"FS", Sev::inform) <<
"\"Plus\" phase viscosity mu_p = " <<
mu_p;
605 <<
"Height of the well in energy functional W = " <<
W;
606 MOFEM_LOG(
"FS", Sev::inform) <<
"Characteristic mesh size h = " <<
h;
607 MOFEM_LOG(
"FS", Sev::inform) <<
"Mobility md = " <<
md;
611 MOFEM_LOG(
"FS", Sev::inform) <<
"Average viscosity mu_ave = " <<
mu_ave;
612 MOFEM_LOG(
"FS", Sev::inform) <<
"Difference viscosity mu_diff = " <<
mu_diff;
642 auto set_problem_bit = [&]() {
665#ifdef PYTHON_INIT_SURFACE
666 auto get_py_surface_init = []() {
667 auto py_surf_init = boost::make_shared<SurfacePython>();
668 CHKERR py_surf_init->surfaceInit(
"surface.py");
669 surfacePythonWeakPtr = py_surf_init;
672 auto py_surf_init = get_py_surface_init();
679 auto dm =
simple->getDM();
683 auto reset_bits = [&]() {
687 start_mask[s] =
true;
689 CHKERR bit_mng->lambdaBitRefLevel(
698 auto add_parent_field = [&](
auto fe,
auto op,
auto field) {
703 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
711 auto h_ptr = boost::make_shared<VectorDouble>();
712 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
713 auto g_ptr = boost::make_shared<VectorDouble>();
714 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
716 auto set_generic = [&](
auto fe) {
718 auto &pip = fe->getOpPtrVector();
726 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
729 fe_parent->getOpPtrVector(), {H1});
735 CHKERR add_parent_field(fe, UDO::OPROW,
"H");
740 CHKERR add_parent_field(fe, UDO::OPROW,
"G");
748 auto post_proc = [&](
auto exe_test) {
750 auto post_proc_fe = boost::make_shared<PostProcEleDomain>(
mField);
751 post_proc_fe->exeTestHook = exe_test;
753 CHKERR set_generic(post_proc_fe);
757 post_proc_fe->getOpPtrVector().push_back(
759 new OpPPMap(post_proc_fe->getPostProcMesh(),
760 post_proc_fe->getMapGaussPts(),
762 {{
"H", h_ptr}, {
"G", g_ptr}},
764 {{
"GRAD_H", grad_h_ptr}, {
"GRAD_G", grad_g_ptr}},
775 CHKERR post_proc_fe->writeFile(
"out_init.h5m");
780 auto solve_init = [&](
auto exe_test) {
783 pip_mng->getOpDomainRhsPipeline().clear();
784 pip_mng->getOpDomainLhsPipeline().clear();
786 auto set_domain_rhs = [&](
auto fe) {
789 auto &pip = fe->getOpPtrVector();
791 CHKERR add_parent_field(fe, UDO::OPROW,
"H");
792 pip.push_back(
new OpRhsH<true>(
"H",
nullptr,
nullptr, h_ptr, grad_h_ptr,
794 CHKERR add_parent_field(fe, UDO::OPROW,
"G");
795 pip.push_back(
new OpRhsG<true>(
"G", h_ptr, grad_h_ptr, g_ptr));
799 auto set_domain_lhs = [&](
auto fe) {
803 auto &pip = fe->getOpPtrVector();
805 CHKERR add_parent_field(fe, UDO::OPROW,
"H");
806 CHKERR add_parent_field(fe, UDO::OPCOL,
"H");
809 CHKERR add_parent_field(fe, UDO::OPCOL,
"G");
812 CHKERR add_parent_field(fe, UDO::OPROW,
"G");
815 CHKERR add_parent_field(fe, UDO::OPCOL,
"H");
820 auto create_subdm = [&]() {
821 auto level_ents_ptr = boost::make_shared<Range>();
827 CHKERR DMSetType(subdm,
"DMMOFEM");
833 for (
auto f : {
"H",
"G"}) {
839 if constexpr (
debug) {
843 dm_ents.subset_by_type(MBVERTEX));
844 dm_ents = subtract(dm_ents, dm_ents.subset_by_type(MBVERTEX));
852 auto subdm = create_subdm();
854 pip_mng->getDomainRhsFE().reset();
855 pip_mng->getDomainLhsFE().reset();
858 pip_mng->getDomainLhsFE()->exeTestHook = exe_test;
859 pip_mng->getDomainRhsFE()->exeTestHook = exe_test;
861 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
862 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
865 auto snes = pip_mng->createSNES(subdm);
868 auto set_section_monitor = [&](
auto solver) {
870 CHKERR SNESMonitorSet(solver,
873 (
void *)(snes_ctx_ptr.get()),
nullptr);
878 auto print_section_field = [&]() {
884 CHKERR PetscSectionGetNumFields(section, &num_fields);
885 for (
int f = 0;
f < num_fields; ++
f) {
889 <<
"Field " <<
f <<
" " << std::string(
field_name);
895 CHKERR set_section_monitor(snes);
896 CHKERR print_section_field();
898 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
902 CHKERR SNESSetFromOptions(snes);
903 CHKERR SNESSolve(snes, PETSC_NULL,
D);
905 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
906 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
917 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
928 pip_mng->getOpDomainRhsPipeline().clear();
929 pip_mng->getOpDomainLhsPipeline().clear();
932 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM_X",
934 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM_X",
936 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM_Y",
938 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM_Y",
940 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX",
"U",
942 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX",
"L",
944 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();
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}
1294 CHKERR prb_mng->buildSubProblem(
"SUB_SOLVER", fields, fields,
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_NULL, 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()) {
1460 for (
auto v : ts_solver_vecs) {
1461 MOFEM_LOG(
"FS", Sev::inform) <<
"Solve projection vector";
1467 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
1468 MOFEM_LOG(
"WORLD", Sev::verbose) <<
"Zero field " << f;
1469 CHKERR field_blas->fieldLambdaOnEntities(zero_dofs, f);
1477 MOFEM_LOG(
"FS", Sev::inform) <<
"Norm V " << get_norm(
v);
1480 CHKERR DMRestoreNamedGlobalVector(dummy_dm,
"TSTheta_X0",
1481 &ts_solver_vecs[0]);
1482 CHKERR DMRestoreNamedGlobalVector(dummy_dm,
"TSTheta_Xdot",
1483 &ts_solver_vecs[1]);
1486 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
1487 MOFEM_LOG(
"WORLD", Sev::verbose) <<
"Zero field " << f;
1488 CHKERR field_blas->fieldLambdaOnEntities(zero_dofs, f);
1497 auto post_proc = [&](
auto exe_test) {
1499 auto post_proc_fe = boost::make_shared<PostProcEleDomain>(
mField);
1500 auto &pip = post_proc_fe->getOpPtrVector();
1501 post_proc_fe->exeTestHook = exe_test;
1509 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1512 fe_parent->getOpPtrVector(), {H1});
1520 auto u_ptr = boost::make_shared<MatrixDouble>();
1521 auto p_ptr = boost::make_shared<VectorDouble>();
1522 auto h_ptr = boost::make_shared<VectorDouble>();
1523 auto g_ptr = boost::make_shared<VectorDouble>();
1525 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"U",
1528 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"P",
1531 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"H",
1534 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"G",
1538 post_proc_fe->getOpPtrVector().push_back(
1540 new OpPPMap(post_proc_fe->getPostProcMesh(),
1541 post_proc_fe->getMapGaussPts(),
1543 {{
"P", p_ptr}, {
"H", h_ptr}, {
"G", g_ptr}},
1555 auto dm =
simple->getDM();
1557 CHKERR post_proc_fe->writeFile(
"out_projection.h5m");
1566 if constexpr (
debug) {
1572 fe_domain_lhs.swap(pip_mng->getDomainLhsFE());
1573 fe_domain_rhs.swap(pip_mng->getDomainRhsFE());
1574 fe_bdy_lhs.swap(pip_mng->getBoundaryLhsFE());
1575 fe_bdy_rhs.swap(pip_mng->getBoundaryRhsFE());
1588 auto add_parent_field_domain = [&](
auto fe,
auto op,
auto field) {
1593 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1601 auto add_parent_field_bdy = [&](
auto fe,
auto op,
auto field) {
1606 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1619 auto dot_u_ptr = boost::make_shared<MatrixDouble>();
1620 auto u_ptr = boost::make_shared<MatrixDouble>();
1621 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
1622 auto dot_h_ptr = boost::make_shared<VectorDouble>();
1623 auto h_ptr = boost::make_shared<VectorDouble>();
1624 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
1625 auto g_ptr = boost::make_shared<VectorDouble>();
1626 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
1627 auto lambda_ptr = boost::make_shared<VectorDouble>();
1628 auto p_ptr = boost::make_shared<VectorDouble>();
1629 auto div_u_ptr = boost::make_shared<VectorDouble>();
1633 auto set_domain_general = [&](
auto fe) {
1635 auto &pip = fe->getOpPtrVector();
1643 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1646 fe_parent->getOpPtrVector(), {H1});
1652 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"U");
1672 "Coordinate system not implemented");
1675 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"H");
1681 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"G");
1686 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"P");
1691 auto set_domain_rhs = [&](
auto fe) {
1693 auto &pip = fe->getOpPtrVector();
1695 CHKERR set_domain_general(fe);
1697 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"U");
1698 pip.push_back(
new OpRhsU(
"U", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr,
1699 grad_h_ptr, g_ptr, p_ptr));
1701 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"H");
1702 pip.push_back(
new OpRhsH<false>(
"H", u_ptr, dot_h_ptr, h_ptr, grad_h_ptr,
1705 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"G");
1706 pip.push_back(
new OpRhsG<false>(
"G", h_ptr, grad_h_ptr, g_ptr));
1708 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"P");
1710 "P", div_u_ptr, [](
const double r,
const double,
const double) {
1714 "P", p_ptr, [](
const double r,
const double,
const double) {
1720 auto set_domain_lhs = [&](
auto fe) {
1722 auto &pip = fe->getOpPtrVector();
1724 CHKERR set_domain_general(fe);
1726 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"U");
1728 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"U");
1729 pip.push_back(
new OpLhsU_dU(
"U", u_ptr, grad_u_ptr, h_ptr));
1730 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"H");
1732 new OpLhsU_dH(
"U",
"H", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr, g_ptr));
1734 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"G");
1735 pip.push_back(
new OpLhsU_dG(
"U",
"G", grad_h_ptr));
1738 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"H");
1740 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"U");
1741 pip.push_back(
new OpLhsH_dU(
"H",
"U", grad_h_ptr));
1742 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"H");
1744 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"G");
1748 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"G");
1750 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"H");
1752 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"G");
1756 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"P");
1758 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"U");
1764 [](
const double r,
const double,
const double) {
1772 [](
const double r,
const double,
const double) {
1779 "Coordinate system not implemented");
1782 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"P");
1783 pip.push_back(
new OpDomainMassP(
"P",
"P", [](
double r,
double,
double) {
1791 auto get_block_name = [](
auto name) {
1792 return boost::format(
"%s(.*)") %
"WETTING_ANGLE";
1795 auto get_blocks = [&](
auto &&name) {
1797 std::regex(name.str()));
1800 auto set_boundary_rhs = [&](
auto fe) {
1802 auto &pip = fe->getOpPtrVector();
1808 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1815 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"U");
1818 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"L");
1823 pip,
mField,
"L", {},
"INFLUX",
1826 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"U");
1829 auto wetting_block = get_blocks(get_block_name(
"WETTING_ANGLE"));
1830 if (wetting_block.size()) {
1837 op_bdy_side->getOpPtrVector(), {H1});
1840 op_bdy_side->getSideFEPtr(),
"", UDO::OPSPACE,
1844 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1847 fe_parent->getOpPtrVector(), {H1});
1853 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPROW,
1855 op_bdy_side->getOpPtrVector().push_back(
1858 pip.push_back(op_bdy_side);
1861 for (
auto &b : wetting_block) {
1863 std::vector<double> attr_vec;
1864 CHKERR b->getMeshsetIdEntitiesByDimension(
1866 b->getAttributes(attr_vec);
1867 if (attr_vec.size() != 1)
1869 "Should be one attribute");
1870 MOFEM_LOG(
"FS", Sev::inform) <<
"Wetting angle: " << attr_vec.front();
1872 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"G");
1874 "G", grad_h_ptr, boost::make_shared<Range>(force_edges),
1882 auto set_boundary_lhs = [&](
auto fe) {
1884 auto &pip = fe->getOpPtrVector();
1890 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1897 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"L");
1898 CHKERR add_parent_field_bdy(fe, UDO::OPCOL,
"U");
1901 auto wetting_block = get_blocks(get_block_name(
"WETTING_ANGLE"));
1902 if (wetting_block.size()) {
1903 auto col_ind_ptr = boost::make_shared<std::vector<VectorInt>>();
1904 auto col_diff_base_ptr = boost::make_shared<std::vector<MatrixDouble>>();
1912 op_bdy_side->getOpPtrVector(), {H1});
1915 op_bdy_side->getSideFEPtr(),
"", UDO::OPSPACE,
1919 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1922 fe_parent->getOpPtrVector(), {H1});
1928 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPROW,
1930 op_bdy_side->getOpPtrVector().push_back(
1932 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPCOL,
1934 op_bdy_side->getOpPtrVector().push_back(
1938 pip.push_back(op_bdy_side);
1941 for (
auto &b : wetting_block) {
1943 std::vector<double> attr_vec;
1944 CHKERR b->getMeshsetIdEntitiesByDimension(
1946 b->getAttributes(attr_vec);
1947 if (attr_vec.size() != 1)
1949 "Should be one attribute");
1951 <<
"wetting angle edges size " << force_edges.size();
1953 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"G");
1955 "G", grad_h_ptr, col_ind_ptr, col_diff_base_ptr,
1956 boost::make_shared<Range>(force_edges), attr_vec.front()));
1965 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
1966 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
1967 CHKERR set_boundary_rhs(pip_mng->getCastBoundaryRhsFE());
1968 CHKERR set_boundary_lhs(pip_mng->getCastBoundaryLhsFE());
1992 boost::shared_ptr<PostProcEleDomainCont> post_proc,
1993 boost::shared_ptr<PostProcEleBdyCont> post_proc_edge,
1994 std::pair<boost::shared_ptr<BoundaryEle>, boost::shared_ptr<VectorDouble>>
2001 MOFEM_LOG(
"FS", Sev::verbose) <<
"Monitor";
2004 MOFEM_LOG(
"FS", Sev::verbose) <<
"Mesh pre proc";
2007 auto post_proc_begin =
2008 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
2010 auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
2018 CHKERR post_proc_end->writeFile(
2019 "out_step_" + boost::lexical_cast<std::string>(
ts_step) +
".h5m");
2020 MOFEM_LOG(
"FS", Sev::verbose) <<
"Mesh pre proc done";
2026 MPI_Allreduce(MPI_IN_PLACE, &(*
liftVec)[0],
SPACE_DIM, MPI_DOUBLE, MPI_SUM,
2030 <<
" lift vec x: " << (*liftVec)[0] <<
" y: " << (*liftVec)[1];
2056 auto setup_subdm = [&](
auto dm) {
2060 auto dm =
simple->getDM();
2061 auto level_ents_ptr = boost::make_shared<Range>();
2062 CHKERR bit_mng->getEntitiesByRefLevel(
2065 CHKERR DMSetType(subdm,
"DMMOFEM");
2071 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
2084 auto get_fe_post_proc = [&](
auto post_proc_mesh) {
2085 auto add_parent_field_domain = [&](
auto fe,
auto op,
auto field) {
2090 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2099 boost::make_shared<PostProcEleDomainCont>(
mField, post_proc_mesh);
2100 post_proc_fe->exeTestHook = [](
FEMethod *fe_ptr) {
2105 auto u_ptr = boost::make_shared<MatrixDouble>();
2106 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
2107 auto h_ptr = boost::make_shared<VectorDouble>();
2108 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
2109 auto p_ptr = boost::make_shared<VectorDouble>();
2110 auto g_ptr = boost::make_shared<VectorDouble>();
2111 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
2114 post_proc_fe->getOpPtrVector(), {H1});
2120 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2123 fe_parent->getOpPtrVector(), {H1});
2129 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"U");
2130 post_proc_fe->getOpPtrVector().push_back(
2132 post_proc_fe->getOpPtrVector().push_back(
2136 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"H");
2137 post_proc_fe->getOpPtrVector().push_back(
2139 post_proc_fe->getOpPtrVector().push_back(
2142 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"P");
2143 post_proc_fe->getOpPtrVector().push_back(
2146 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"G");
2147 post_proc_fe->getOpPtrVector().push_back(
2149 post_proc_fe->getOpPtrVector().push_back(
2154 post_proc_fe->getOpPtrVector().push_back(
2157 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
2159 {{
"H", h_ptr}, {
"P", p_ptr}, {
"G", g_ptr}},
2161 {{
"U", u_ptr}, {
"H_GRAD", grad_h_ptr}, {
"G_GRAD", grad_g_ptr}},
2163 {{
"GRAD_U", grad_u_ptr}},
2171 return post_proc_fe;
2174 auto get_bdy_post_proc_fe = [&](
auto post_proc_mesh) {
2175 auto add_parent_field_bdy = [&](
auto fe,
auto op,
auto field) {
2176 return setParentDofs(
2180 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2189 boost::make_shared<PostProcEleBdyCont>(mField, post_proc_mesh);
2190 post_proc_fe->exeTestHook = [](
FEMethod *fe_ptr) {
2199 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2208 OpGetNormal(boost::shared_ptr<VectorDouble> l_ptr,
2209 boost::shared_ptr<MatrixDouble> n_ptr)
2219 auto t_n = getFTensor1FromMat<SPACE_DIM>(*ptrNormal);
2220 for (
auto gg = 0; gg !=
getGaussPts().size2(); ++gg) {
2221 t_n(
i) = t_n_fe(
i) * t_l / std::sqrt(t_n_fe(
i) * t_n_fe(
i));
2230 boost::shared_ptr<VectorDouble> ptrL;
2231 boost::shared_ptr<MatrixDouble> ptrNormal;
2234 auto u_ptr = boost::make_shared<MatrixDouble>();
2235 auto p_ptr = boost::make_shared<VectorDouble>();
2236 auto lambda_ptr = boost::make_shared<VectorDouble>();
2237 auto normal_l_ptr = boost::make_shared<MatrixDouble>();
2239 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW,
"U");
2240 post_proc_fe->getOpPtrVector().push_back(
2243 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW,
"L");
2244 post_proc_fe->getOpPtrVector().push_back(
2246 post_proc_fe->getOpPtrVector().push_back(
2247 new OpGetNormal(lambda_ptr, normal_l_ptr));
2249 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW,
"P");
2250 post_proc_fe->getOpPtrVector().push_back(
2254 op_ptr->doWorkRhsHook = [&](
DataOperator *base_op_ptr,
int side,
2264 post_proc_fe->getOpPtrVector().push_back(
2266 new OpPPMap(post_proc_fe->getPostProcMesh(),
2267 post_proc_fe->getMapGaussPts(),
2281 return post_proc_fe;
2284 auto get_lift_fe = [&]() {
2285 auto add_parent_field_bdy = [&](
auto fe,
auto op,
auto field) {
2286 return setParentDofs(
2290 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2298 auto fe = boost::make_shared<BoundaryEle>(mField);
2299 fe->exeTestHook = [](
FEMethod *fe_ptr) {
2304 auto lift_ptr = boost::make_shared<VectorDouble>();
2305 auto p_ptr = boost::make_shared<VectorDouble>();
2306 auto ents_ptr = boost::make_shared<Range>();
2312 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2319 std::vector<const CubitMeshSets *> vec_ptr;
2321 std::regex(
"LIFT"), vec_ptr);
2322 for (
auto m_ptr : vec_ptr) {
2327 ents_ptr->merge(ents);
2330 MOFEM_LOG(
"FS", Sev::noisy) <<
"Lift ents " << (*ents_ptr);
2332 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"P");
2333 fe->getOpPtrVector().push_back(
2335 fe->getOpPtrVector().push_back(
2338 return std::make_pair(fe, lift_ptr);
2341 auto set_post_proc_monitor = [&](
auto ts) {
2345 boost::shared_ptr<FEMethod> null_fe;
2346 auto post_proc_mesh = boost::make_shared<moab::Core>();
2347 auto monitor_ptr = boost::make_shared<Monitor>(
2349 get_fe_post_proc(post_proc_mesh), get_bdy_post_proc_fe(post_proc_mesh),
2352 null_fe, monitor_ptr);
2356 auto dm =
simple->getDM();
2360 auto ts_pre_post_proc = boost::make_shared<TSPrePostProc>();
2365 ptr->fsRawPtr =
this;
2366 ptr->solverSubDM = create_solver_dm(
simple->getDM());
2370 CHKERR VecAssemblyBegin(ptr->globSol);
2371 CHKERR VecAssemblyEnd(ptr->globSol);
2373 auto sub_ts = pip_mng->createTSIM(ptr->solverSubDM);
2375 CHKERR set_post_proc_monitor(sub_ts);
2379 CHKERR TSSetFromOptions(ts);
2383 auto print_fields_in_section = [&]() {
2386 simple->getProblemName());
2387 PetscInt num_fields;
2388 CHKERR PetscSectionGetNumFields(section, &num_fields);
2389 for (
int f = 0;
f < num_fields; ++
f) {
2393 <<
"Field " <<
f <<
" " << std::string(
field_name);
2398 CHKERR print_fields_in_section();
2400 CHKERR TSSolve(ts, ptr->globSol);
2408#ifdef PYTHON_INIT_SURFACE
2413 const char param_file[] =
"param_file.petsc";
2417 auto core_log = logging::core::get();
2425 DMType dm_name =
"DMMOFEM";
2430 moab::Core mb_instance;
2431 moab::Interface &moab = mb_instance;
2448#ifdef PYTHON_INIT_SURFACE
2449 if (Py_FinalizeEx() < 0) {
2465 "can not get vertices on bit 0");
2470 Range plus_range, minus_range;
2471 std::vector<EntityHandle> plus, minus;
2474 for (
auto p = vertices.pair_begin();
p != vertices.pair_end(); ++
p) {
2476 const auto f =
p->first;
2477 const auto s =
p->second;
2482 auto it = dofs_mi.lower_bound(lo_uid);
2483 const auto hi_it = dofs_mi.upper_bound(hi_uid);
2487 plus.reserve(std::distance(it, hi_it));
2488 minus.reserve(std::distance(it, hi_it));
2490 for (; it != hi_it; ++it) {
2491 const auto v = (*it)->getFieldData();
2493 plus.push_back((*it)->getEnt());
2495 minus.push_back((*it)->getEnt());
2498 plus_range.insert_list(plus.begin(), plus.end());
2499 minus_range.insert_list(minus.begin(), minus.end());
2504 <<
"Plus range " << plus_range << endl;
2506 <<
"Minus range " << minus_range << endl;
2509 auto get_elems = [&](
auto &ents,
auto bit,
auto mask) {
2512 moab::Interface::UNION),
2513 "can not get adjacencies");
2515 "can not filter elements with bit 0");
2519 CHKERR comm_mng->synchroniseEntities(plus_range);
2520 CHKERR comm_mng->synchroniseEntities(minus_range);
2523 ele_plus[0] = get_elems(plus_range,
bit(0),
BitRefLevel().set());
2524 ele_minus[0] = get_elems(minus_range,
bit(0),
BitRefLevel().set());
2525 auto common = intersect(ele_plus[0], ele_minus[0]);
2526 ele_plus[0] = subtract(ele_plus[0], common);
2527 ele_minus[0] = subtract(ele_minus[0], common);
2529 auto get_children = [&](
auto &
p,
auto &
c) {
2531 CHKERR bit_mng->updateRangeByChildren(
p,
c);
2543 auto get_level = [&](
auto &
p,
auto &
m,
auto z,
auto bit,
auto mask) {
2546 bit_mng->getEntitiesByDimAndRefLevel(
bit, mask,
SPACE_DIM,
l),
2547 "can not get vertices on bit");
2550 for (
auto f = 0; f != z; ++f) {
2554 l = get_elems(conn,
bit, mask);
2559 std::vector<Range> vec_levels(
nb_levels);
2561 vec_levels[
l] = get_level(ele_plus[
l], ele_minus[
l], 2 * overlap,
bit(
l),
2565 if constexpr (
debug) {
2568 std::string name = (boost::format(
"out_r%d.h5m") %
l).str();
2587 start_mask[s] =
true;
2593 Range prev_level_skin;
2597 CHKERR bit_mng->lambdaBitRefLevel(
2604 auto set_levels = [&](
auto &&
2614 auto get_adj = [&](
auto ents) {
2620 ents.subset_by_dimension(
SPACE_DIM), d,
false, conn,
2621 moab::Interface::UNION),
2635 CHKERR bit_mng->updateRangeByParent(vec_levels[
l], parents);
2638 level_prev = subtract(level_prev, parents);
2640 level_prev.merge(vec_levels[
l]);
2648 CHKERR bit_mng->getEntitiesByDimAndRefLevel(
2650 level = get_adj(level);
2659 auto set_skins = [&]() {
2663 ParallelComm *pcomm =
2672 "can't get bit level");
2679 auto get_level_skin = [&]() {
2685 CHKERR pcomm->filter_pstatus(skin_level_mesh,
2686 PSTATUS_SHARED | PSTATUS_MULTISHARED,
2687 PSTATUS_NOT, -1,
nullptr);
2690 skin_level = subtract(skin_level,
2694 Range skin_level_verts;
2697 skin_level.merge(skin_level_verts);
2700 bit_prev.set(
l - 1);
2704 skin.merge(subtract(skin_level, level_prev));
2710 auto resolve_shared = [&](
auto &&skin) {
2711 Range tmp_skin = skin;
2713 map<int, Range> map_procs;
2715 tmp_skin, &map_procs);
2717 Range from_other_procs;
2718 for (
auto &
m : map_procs) {
2720 from_other_procs.merge(
m.second);
2724 auto common = intersect(
2725 skin, from_other_procs);
2726 skin.merge(from_other_procs);
2730 if (!common.empty()) {
2732 skin = subtract(skin, common);
2739 auto get_parent_level_skin = [&](
auto skin) {
2741 CHKERR bit_mng->updateRangeByParent(
2742 skin.subset_by_dimension(
SPACE_DIM - 1), skin_parents);
2743 Range skin_parent_verts;
2746 skin_parents.merge(skin_parent_verts);
2749 return skin_parents;
2752 auto child_skin = resolve_shared(get_level_skin());
2753 auto parent_skin = get_parent_level_skin(child_skin);
2755 child_skin = subtract(child_skin, parent_skin);
2763 auto set_current = [&]() {
2772 last_level = subtract(last_level, skin_child);
2784 if constexpr (
debug) {
2787 std::string name = (boost::format(
"out_level%d.h5m") %
l).str();
2790 "PARALLEL=WRITE_PART");
2794 "MOAB",
"PARALLEL=WRITE_PART");
2797 "MOAB",
"PARALLEL=WRITE_PART");
2801 "MOAB",
"PARALLEL=WRITE_PART");
2804 "MOAB",
"PARALLEL=WRITE_PART");
2812 boost::shared_ptr<FEMethod> fe_top, std::string
field_name,
2813 ForcesAndSourcesCore::UserDataOperator::OpType op,
2815 boost::function<boost::shared_ptr<ForcesAndSourcesCore>()> get_elem,
2822 boost::function<void(boost::shared_ptr<ForcesAndSourcesCore>,
int)>
2824 [&](boost::shared_ptr<ForcesAndSourcesCore> parent_fe_pt,
int level) {
2829 auto fe_ptr_current = get_elem();
2833 boost::dynamic_pointer_cast<ForcesAndSourcesCore>(
2838 if (op == ForcesAndSourcesCore::UserDataOperator::OPSPACE) {
2841 parent_fe_pt->getOpPtrVector().push_back(
2845 H1, op, fe_ptr_current,
2856 parent_fe_pt->getOpPtrVector().push_back(
2871 add_parent_level(boost::dynamic_pointer_cast<ForcesAndSourcesCore>(fe_top),
2886 auto cut_off_dofs = [&]() {
2889 auto &m_field = ptr->fsRawPtr->mField;
2891 Range current_verts;
2893 CHKERR bit_mng->getEntitiesByTypeAndRefLevel(
2896 auto cut_off_verts = [&](boost::shared_ptr<FieldEntity> ent_ptr) {
2898 for (
auto &
h : ent_ptr->getEntFieldData()) {
2904 auto field_blas = m_field.getInterface<
FieldBlas>();
2905 CHKERR field_blas->fieldLambdaOnEntities(cut_off_verts,
"H",
2914 MOFEM_LOG(
"FS", Sev::inform) <<
"Run step pre proc";
2916 auto &m_field = ptr->fsRawPtr->mField;
2926 auto get_norm = [&](
auto x) {
2928 CHKERR VecNorm(x, NORM_2, &nrm);
2933 auto refine_problem = [&]() {
2935 MOFEM_LOG(
"FS", Sev::inform) <<
"Refine problem";
2937 CHKERR ptr->fsRawPtr->projectData();
2943 auto set_jacobian_operators = [&]() {
2946 CHKERR KSPReset(ptr->subKSP);
2951 auto set_solution = [&]() {
2953 MOFEM_LOG(
"FS", Sev::inform) <<
"Set solution";
2955 PetscObjectState state;
2966 INSERT_VALUES, SCATTER_FORWARD);
2969 <<
"Set solution, vector norm " << get_norm(ptr->globSol);
2974 PetscObjectTypeCompare((PetscObject)ts, TSTHETA, &is_theta);
2978 CHKERR set_jacobian_operators();
2983 "Sorry, only TSTheta handling is implemented");
2988 PetscBarrier((PetscObject)ts);
3000 auto &m_field = ptr->fsRawPtr->mField;
3012 auto sub_u = ptr->getSubVector();
3015 auto scatter = ptr->getScatter(sub_u, u,
R);
3017 auto apply_scatter_and_update = [&](
auto x,
auto sub_x) {
3019 CHKERR VecScatterBegin(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3020 CHKERR VecScatterEnd(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3021 CHKERR VecGhostUpdateBegin(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3022 CHKERR VecGhostUpdateEnd(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3026 CHKERR apply_scatter_and_update(u, sub_u);
3027 CHKERR apply_scatter_and_update(u_t, sub_u_t);
3030 CHKERR VecScatterBegin(scatter, sub_f, f, INSERT_VALUES, SCATTER_FORWARD);
3031 CHKERR VecScatterEnd(scatter, sub_f, f, INSERT_VALUES, SCATTER_FORWARD);
3037 PetscReal a, Mat
A, Mat
B,
3041 auto sub_u = ptr->getSubVector();
3043 auto scatter = ptr->getScatter(sub_u, u,
R);
3045 auto apply_scatter_and_update = [&](
auto x,
auto sub_x) {
3047 CHKERR VecScatterBegin(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3048 CHKERR VecScatterEnd(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3049 CHKERR VecGhostUpdateBegin(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3050 CHKERR VecGhostUpdateEnd(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3054 CHKERR apply_scatter_and_update(u, sub_u);
3055 CHKERR apply_scatter_and_update(u_t, sub_u_t);
3058 ptr->tsCtxPtr.get());
3067 auto get_norm = [&](
auto x) {
3069 CHKERR VecNorm(x, NORM_2, &nrm);
3073 auto sub_u = ptr->getSubVector();
3074 auto scatter = ptr->getScatter(sub_u, u,
R);
3075 CHKERR VecScatterBegin(scatter, u, sub_u, INSERT_VALUES, SCATTER_REVERSE);
3076 CHKERR VecScatterEnd(scatter, u, sub_u, INSERT_VALUES, SCATTER_REVERSE);
3077 CHKERR VecGhostUpdateBegin(sub_u, INSERT_VALUES, SCATTER_FORWARD);
3078 CHKERR VecGhostUpdateEnd(sub_u, INSERT_VALUES, SCATTER_FORWARD);
3081 <<
"u norm " << get_norm(u) <<
" u sub nom " << get_norm(sub_u);
3091 MOFEM_LOG(
"FS", Sev::verbose) <<
"SetUP sub PC";
3092 ptr->subKSP =
createKSP(ptr->fsRawPtr->mField.get_comm());
3093 CHKERR KSPSetFromOptions(ptr->subKSP);
3101 auto sub_x = ptr->getSubVector();
3103 auto scatter = ptr->getScatter(sub_x, pc_x,
R);
3104 CHKERR VecScatterBegin(scatter, pc_f, sub_f, INSERT_VALUES,
3106 CHKERR VecScatterEnd(scatter, pc_f, sub_f, INSERT_VALUES, SCATTER_REVERSE);
3107 CHKERR KSPSetOperators(ptr->subKSP, ptr->subB, ptr->subB);
3108 MOFEM_LOG(
"FS", Sev::verbose) <<
"PCShell solve";
3109 CHKERR KSPSolve(ptr->subKSP, sub_f, sub_x);
3110 MOFEM_LOG(
"FS", Sev::verbose) <<
"PCShell solve <- done";
3111 CHKERR VecScatterBegin(scatter, sub_x, pc_x, INSERT_VALUES,
3113 CHKERR VecScatterEnd(scatter, sub_x, pc_x, INSERT_VALUES, SCATTER_FORWARD);
3120 auto prb_ptr = ptr->fsRawPtr->mField.get_problem(
"SUB_SOLVER");
3121 if (
auto sub_data = prb_ptr->getSubData()) {
3122 auto is = sub_data->getSmartColIs();
3149 auto dm =
simple->getDM();
3157 CHKERR TSGetSNES(ts, &snes);
3160 auto set_section_monitor = [&](
auto snes) {
3162 CHKERR SNESMonitorSet(snes,
3165 (
void *)(snes_ctx_ptr.get()),
nullptr);
3169 CHKERR set_section_monitor(snes);
3171 auto ksp =
createKSP(m_field.get_comm());
3172 CHKERR KSPSetType(ksp, KSPPREONLY);
3174 CHKERR PCSetType(sub_pc, PCSHELL);
3177 CHKERR KSPSetPC(ksp, sub_pc);
3178 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 >::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
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
constexpr int U_FIELD_DIM
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
static boost::weak_ptr< TSPrePostProc > tsPrePostProc
FTensor::Index< 'j', SPACE_DIM > j
double lambda
surface tension
FTensor::Index< 'k', SPACE_DIM > k
OpDomainMassH OpDomainMassG
FTensor::Index< 'i', SPACE_DIM > i
auto init_h
Initialisation function.
FTensor::Index< 'l', SPACE_DIM > l
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
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpDomainAssembleScalar
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
auto wetting_angle_sub_stepping
FormsIntegrators< BoundaryEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, 1 > OpBoundaryMassL
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
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
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
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)
PipelineManager::ElementsAndOpsByDim< SPACE_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 jacobina 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 *snes_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.
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
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
virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
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
default operator for TRI element
auto getFTensor1NormalsAtGaussPts()
get normal at integration points
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
@ OPSPACE
operator do Work is execute on space data
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
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 filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Approximate field values for given petsc vector.
Get values at integration pts for tensor filed 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< VectorDouble > > DataMapVec
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
Operator to execute finite element instance on parent element. This operator is typically used to pro...
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
PipelineManager interface.
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
intrusive_ptr for managing petsc objects
PetscInt ts_step
time step number
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce 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