Using PipelineManager interface calculate the divergence of base functions, and integral of flux on the boundary. Since the h-div space is used, volume integral and boundary integral should give the same result.
#include <petsc/private/petscimpl.h>
static char help[] =
"...\n\n";
#ifdef PYTHON_INIT_SURFACE
#include <boost/python.hpp>
#include <boost/python/def.hpp>
namespace bp = boost::python;
struct SurfacePython {
SurfacePython() = default;
virtual ~SurfacePython() = default;
try {
auto main_module = bp::import("__main__");
mainNamespace = main_module.attr("__dict__");
bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
surfaceFun = mainNamespace["surface"];
} catch (bp::error_already_set const &) {
PyErr_Print();
}
};
double x,
double y,
double z,
double eta,
double &s
) {
try {
s = bp::extract<double>(surfaceFun(x, y, z,
eta));
} catch (bp::error_already_set const &) {
PyErr_Print();
}
}
private:
bp::object mainNamespace;
bp::object surfaceFun;
};
static boost::weak_ptr<SurfacePython> surfacePythonWeakPtr;
#endif
template <int DIM>
A>::LinearForm<I>::OpBaseTimesScalar<
BASE_DIM>;
A>::LinearForm<I>::OpBaseTimesScalar<
BASE_DIM>;
template <CoordinateTypes COORD_TYPE>
constexpr
bool debug =
true;
};
};
double mu_m = 0.010101 * 1e2;
double mu_p = 0.000182 * 1e2;
double tol = std::numeric_limits<float>::epsilon();
else
return 1.;
};
constexpr int sub_stepping = 16;
return std::min(1., static_cast<double>(ts_step) / sub_stepping);
};
auto my_max = [](
const double x) {
return (x - 1 + std::abs(x + 1)) / 2; };
auto my_min = [](
const double x) {
return (x + 1 - std::abs(x - 1)) / 2; };
return 1.;
else
return 0.;
};
auto phase_function = [](
const double h,
const double diff,
const double ave) {
};
return diff;
};
auto get_f = [](
const double h) {
return 4 *
W *
h * (
h *
h - 1); };
auto get_f_dh = [](
const double h) {
return 4 *
W * (3 *
h *
h - 1); };
};
const double h3 = h2 *
h;
return md * (2 * h3 - 3 * h2 + 1);
else
return md * (-2 * h3 - 3 * h2 + 1);
};
return md * (6 *
h * (
h - 1));
else
return md * (-6 *
h * (
h + 1));
};
auto get_D = [](
const double A) {
return t_D;
};
constexpr
double R = 0.0125;
constexpr
double A =
R * 0.2;
const double theta = atan2(
r, y);
const double w =
R +
A * cos(
n * theta);
const double d = std::sqrt(
r *
r + y * y);
return tanh((
w -
d) / (
eta * std::sqrt(2)));
};
constexpr double y0 = 0.5;
constexpr
double R = 0.4;
y -= y0;
const double d = std::sqrt(
r *
r + y * y);
return tanh((
R -
d) / (
eta * std::sqrt(2)));
};
constexpr double water_height = 0.;
return tanh((water_height - y) / (
eta * std::sqrt(2)));
;
};
return -tanh((-0.039 - x) / (
eta * std::sqrt(2)));
};
auto init_h = [](
double r,
double y,
double theta) {
#ifdef PYTHON_INIT_SURFACE
double s = 1;
if (auto ptr = surfacePythonWeakPtr.lock()) {
"error eval python");
}
return s;
#else
#endif
};
auto wetting_angle = [](
double water_level) {
return water_level; };
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel();
};
int g_size;
MPI_Allreduce(&l_size, &g_size, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
return g_size;
};
CHKERR moab.add_entities(*out_meshset,
r);
CHKERR moab.write_file(name.c_str(),
"MOAB",
"PARALLEL=WRITE_PART",
out_meshset->get_ptr(), 1);
}
};
std::vector<EntityHandle> ents_vec;
auto dofs = prb_ptr->numeredRowDofsPtr;
ents_vec.reserve(std::distance(lo_it, hi_it));
for (; lo_it != hi_it; ++lo_it) {
ents_vec.push_back((*lo_it)->getEnt());
}
std::sort(ents_vec.begin(), ents_vec.end());
auto it = std::unique(ents_vec.begin(), ents_vec.end());
r.insert_list(ents_vec.begin(), it);
};
std::vector<EntityHandle> ents_vec;
auto dofs = prb_ptr->numeredRowDofsPtr;
ents_vec.reserve(dofs->size());
ents_vec.push_back(
d->getEnt());
}
std::sort(ents_vec.begin(), ents_vec.end());
auto it = std::unique(ents_vec.begin(), ents_vec.end());
r.insert_list(ents_vec.begin(), it);
};
private:
void *ctx);
PetscReal
a, Mat
A, Mat B,
void *ctx);
void *ctx);
boost::shared_ptr<SnesCtx>
snesCtxPtr;
boost::shared_ptr<TsCtx>
tsCtxPtr;
};
private:
std::vector<Range> findEntitiesCrossedByPhaseInterface(size_t overlap);
boost::shared_ptr<FEMethod> fe_top, std::string
field_name,
boost::function<boost::shared_ptr<ForcesAndSourcesCore>()> get_elem,
};
}
MOFEM_LOG(
"FS", Sev::inform) <<
"Read mesh for problem";
simple->getParentAdjacencies() =
true;
}
PETSC_NULL);
PETSC_NULL);
const char *coord_type_names[] = {"cartesian", "polar", "cylindrical",
"spherical"};
<<
"Number of refinement levels nb_levels = " <<
nb_levels;
"-rho_m", &
rho_m, PETSC_NULL);
"Height of the well in energy functional", "-W",
MOFEM_LOG(
"FS", Sev::inform) <<
"Acceleration a0 = " <<
a0;
MOFEM_LOG(
"FS", Sev::inform) <<
"\"Minus\" phase density rho_m = " <<
rho_m;
MOFEM_LOG(
"FS", Sev::inform) <<
"\"Minus\" phase viscosity mu_m = " <<
mu_m;
MOFEM_LOG(
"FS", Sev::inform) <<
"\"Plus\" phase density rho_p = " <<
rho_p;
MOFEM_LOG(
"FS", Sev::inform) <<
"\"Plus\" phase viscosity mu_p = " <<
mu_p;
<<
"Height of the well in energy functional W = " <<
W;
MOFEM_LOG(
"FS", Sev::inform) <<
"Characteristic mesh size h = " <<
h;
MOFEM_LOG(
"FS", Sev::inform) <<
"Mobility md = " <<
md;
auto set_problem_bit = [&]() {
};
}
#ifdef PYTHON_INIT_SURFACE
auto get_py_surface_init = []() {
auto py_surf_init = boost::make_shared<SurfacePython>();
CHKERR py_surf_init->surfaceInit(
"surface.py");
surfacePythonWeakPtr = py_surf_init;
return py_surf_init;
};
auto py_surf_init = get_py_surface_init();
#endif
auto reset_bits = [&]() {
start_mask[s] = true;
CHKERR bit_mng->lambdaBitRefLevel(
};
auto add_parent_field = [&](auto fe, auto op, auto field) {
return setParentDofs(
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
};
auto h_ptr = boost::make_shared<VectorDouble>();
auto grad_h_ptr = boost::make_shared<MatrixDouble>();
auto g_ptr = boost::make_shared<VectorDouble>();
auto grad_g_ptr = boost::make_shared<MatrixDouble>();
auto set_generic = [&](auto fe) {
auto &pip = fe->getOpPtrVector();
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
fe_parent->getOpPtrVector(), {H1});
return fe_parent;
},
CHKERR add_parent_field(fe, UDO::OPROW,
"H");
pip.push_back(
CHKERR add_parent_field(fe, UDO::OPROW,
"G");
pip.push_back(
};
auto post_proc = [&](auto exe_test) {
auto post_proc_fe = boost::make_shared<PostProcEleDomain>(mField);
post_proc_fe->exeTestHook = exe_test;
CHKERR set_generic(post_proc_fe);
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
{{"H", h_ptr}, {"G", g_ptr}},
{{"GRAD_H", grad_h_ptr}, {"GRAD_G", grad_g_ptr}},
{},
{}
)
);
CHKERR post_proc_fe->writeFile(
"out_init.h5m");
};
auto solve_init = [&](auto exe_test) {
pip_mng->getOpDomainRhsPipeline().clear();
pip_mng->getOpDomainLhsPipeline().clear();
auto set_domain_rhs = [&](auto fe) {
auto &pip = fe->getOpPtrVector();
CHKERR add_parent_field(fe, UDO::OPROW,
"H");
pip.push_back(
new OpRhsH<true>(
"H",
nullptr,
nullptr, h_ptr, grad_h_ptr,
grad_g_ptr));
CHKERR add_parent_field(fe, UDO::OPROW,
"G");
pip.push_back(
new OpRhsG<true>(
"G", h_ptr, grad_h_ptr, g_ptr));
};
auto set_domain_lhs = [&](auto fe) {
auto &pip = fe->getOpPtrVector();
CHKERR add_parent_field(fe, UDO::OPROW,
"H");
CHKERR add_parent_field(fe, UDO::OPCOL,
"H");
CHKERR add_parent_field(fe, UDO::OPCOL,
"G");
CHKERR add_parent_field(fe, UDO::OPROW,
"G");
CHKERR add_parent_field(fe, UDO::OPCOL,
"H");
};
auto create_subdm = [&]() {
auto level_ents_ptr = boost::make_shared<Range>();
DM subdm;
CHKERR DMSetType(subdm,
"DMMOFEM");
for (
auto f : {
"H",
"G"}) {
}
dm_ents.subset_by_type(MBVERTEX));
dm_ents = subtract(dm_ents, dm_ents.subset_by_type(MBVERTEX));
}
}
};
auto subdm = create_subdm();
pip_mng->getDomainRhsFE().reset();
pip_mng->getDomainLhsFE().reset();
pip_mng->getDomainLhsFE()->exeTestHook = exe_test;
pip_mng->getDomainRhsFE()->exeTestHook = exe_test;
CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
auto snes = pip_mng->createSNES(subdm);
auto set_section_monitor = [&](auto solver) {
(void *)(snes_ctx_ptr.get()), nullptr);
};
auto print_section_field = [&]() {
auto section =
PetscInt num_fields;
CHKERR PetscSectionGetNumFields(section, &num_fields);
for (
int f = 0;
f < num_fields; ++
f) {
}
};
CHKERR set_section_monitor(snes);
for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
}
CHKERR SNESSetFromOptions(snes);
CHKERR SNESSolve(snes, PETSC_NULL,
D);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
};
for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
}
});
});
pip_mng->getOpDomainRhsPipeline().clear();
pip_mng->getOpDomainLhsPipeline().clear();
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM_X",
"U", 0, 0);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM_X",
"L", 0, 0);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM_Y",
"U", 1, 1);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM_Y",
"L", 1, 1);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX",
"U",
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX",
"L",
0, 0);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"ZERO",
"L", 0, 0);
}
auto fe_domain_lhs = pip_mng->getDomainLhsFE();
auto fe_domain_rhs = pip_mng->getDomainRhsFE();
auto fe_bdy_lhs = pip_mng->getBoundaryLhsFE();
auto fe_bdy_rhs = pip_mng->getBoundaryRhsFE();
pip_mng->getDomainLhsFE().reset();
pip_mng->getDomainRhsFE().reset();
pip_mng->getBoundaryLhsFE().reset();
pip_mng->getBoundaryRhsFE().reset();
auto add_parent_field_domain = [&](
auto fe,
auto op,
auto field,
auto bit) {
return setParentDofs(
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
};
auto add_parent_field_bdy = [&](
auto fe,
auto op,
auto field,
auto bit) {
return setParentDofs(
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
};
auto create_run_parent_op = [&](auto parent_fe_ptr, auto this_fe_ptr,
auto fe_bit) {
auto parent_mask = fe_bit;
parent_mask.flip();
Sev::inform);
};
auto get_parents_vel_fe_ptr = [&](auto this_fe_ptr, auto fe_bit) {
std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
parents_elems_ptr_vec.emplace_back(
boost::make_shared<DomianParentEle>(mField));
parents_elems_ptr_vec[
l - 1]->getOpPtrVector().push_back(
create_run_parent_op(parents_elems_ptr_vec[
l], this_fe_ptr, fe_bit));
}
return parents_elems_ptr_vec[0];
};
auto solve_projection = [&](auto exe_test) {
auto set_domain_rhs = [&](auto fe) {
auto &pip = fe->getOpPtrVector();
auto u_ptr = boost::make_shared<MatrixDouble>();
auto p_ptr = boost::make_shared<VectorDouble>();
auto h_ptr = boost::make_shared<VectorDouble>();
auto g_ptr = boost::make_shared<VectorDouble>();
auto eval_fe_ptr = boost::make_shared<DomianParentEle>(mField);
{
auto &pip = eval_fe_ptr->getOpPtrVector();
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW,
"U",
CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW,
"P",
CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW,
"H",
CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW,
"G",
}
auto parent_eval_fe_ptr =
pip.push_back(create_run_parent_op(parent_eval_fe_ptr, eval_fe_ptr,
auto assemble_fe_ptr = boost::make_shared<DomianParentEle>(mField);
{
auto &pip = assemble_fe_ptr->getOpPtrVector();
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW,
"U",
CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW,
"P",
CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW,
"H",
CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW,
"G",
}
auto parent_assemble_fe_ptr =
pip.push_back(create_run_parent_op(
};
auto set_domain_lhs = [&](auto fe) {
auto &pip = fe->getOpPtrVector();
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
CHKERR add_parent_field_domain(fe, UDO::OPROW,
"U",
CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"U",
CHKERR add_parent_field_domain(fe, UDO::OPROW,
"P",
CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"P",
CHKERR add_parent_field_domain(fe, UDO::OPROW,
"H",
CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"H",
CHKERR add_parent_field_domain(fe, UDO::OPROW,
"G",
CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"G",
};
auto set_bdy_rhs = [&](auto fe) {
auto &pip = fe->getOpPtrVector();
auto l_ptr = boost::make_shared<VectorDouble>();
auto eval_fe_ptr = boost::make_shared<BoundaryParentEle>(mField);
{
auto &pip = eval_fe_ptr->getOpPtrVector();
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
CHKERR add_parent_field_bdy(eval_fe_ptr, UDO::OPROW,
"L",
}
auto parent_eval_fe_ptr =
pip.push_back(create_run_parent_op(parent_eval_fe_ptr, eval_fe_ptr,
auto assemble_fe_ptr = boost::make_shared<BoundaryParentEle>(mField);
{
auto &pip = assemble_fe_ptr->getOpPtrVector();
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
OpLSize(boost::shared_ptr<VectorDouble> l_ptr)
if (lPtr->size() != getGaussPts().size2()) {
lPtr->resize(getGaussPts().size2());
lPtr->clear();
}
}
private:
boost::shared_ptr<VectorDouble> lPtr;
};
pip.push_back(new OpLSize(l_ptr));
CHKERR add_parent_field_bdy(assemble_fe_ptr, UDO::OPROW,
"L",
}
auto parent_assemble_fe_ptr =
pip.push_back(create_run_parent_op(
};
auto set_bdy_lhs = [&](auto fe) {
auto &pip = fe->getOpPtrVector();
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"L",
CHKERR add_parent_field_bdy(fe, UDO::OPCOL,
"L",
};
auto level_ents_ptr = boost::make_shared<Range>();
auto get_prj_ents = [&]() {
auto common_ents = intersect(prj_mesh, *level_ents_ptr);
prj_mesh = subtract(unite(*level_ents_ptr, prj_mesh), common_ents)
return prj_mesh;
};
auto prj_ents = get_prj_ents();
auto rebuild = [&]() {
std::vector<std::string> fields{"U", "P", "H", "G", "L"};
std::map<std::string, boost::shared_ptr<Range>> range_maps{
{"U", level_ents_ptr},
{"P", level_ents_ptr},
{"H", level_ents_ptr},
{"G", level_ents_ptr},
{"L", level_ents_ptr}
};
CHKERR prb_mng->buildSubProblem(
"SUB_SOLVER", fields, fields,
simple->getProblemName(), PETSC_TRUE,
&range_maps, &range_maps);
CHKERR prb_mng->partitionFiniteElements(
"SUB_SOLVER",
true, 0,
CHKERR prb_mng->partitionGhostDofsOnDistributedMesh(
"SUB_SOLVER");
};
MOFEM_LOG(
"FS", Sev::verbose) <<
"Create projection problem";
DM subdm;
CHKERR DMSetType(subdm,
"DMMOFEM");
}
MOFEM_LOG(
"FS", Sev::inform) <<
"Nothing to project";
};
auto create_dummy_dm = [&]() {
simple->getProblemName().c_str(),
"create dummy dm");
return dummy_dm;
};
auto subdm = create_subdm();
if (subdm) {
pip_mng->getDomainRhsFE().reset();
pip_mng->getDomainLhsFE().reset();
pip_mng->getBoundaryRhsFE().reset();
pip_mng->getBoundaryLhsFE().reset();
pip_mng->getDomainLhsFE()->exeTestHook = exe_test;
pip_mng->getDomainRhsFE()->exeTestHook = [](
FEMethod *fe_ptr) {
};
pip_mng->getBoundaryLhsFE()->exeTestHook = exe_test;
pip_mng->getBoundaryRhsFE()->exeTestHook = [](
FEMethod *fe_ptr) {
};
CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
CHKERR set_bdy_rhs(pip_mng->getCastBoundaryRhsFE());
CHKERR set_bdy_lhs(pip_mng->getCastBoundaryLhsFE());
auto ksp = pip_mng->createKSP(subdm);
CHKERR KSPSetFromOptions(ksp);
auto get_norm = [&](auto x) {
double nrm;
CHKERR VecNorm(x, NORM_2, &nrm);
return nrm;
};
auto zero_dofs = [](boost::shared_ptr<FieldEntity> ent_ptr) {
for (
auto &
v : ent_ptr->getEntFieldData()) {
}
};
auto solve = [&](
auto S) {
CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(S, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(S, INSERT_VALUES, SCATTER_FORWARD);
};
MOFEM_LOG(
"FS", Sev::inform) <<
"Solve projection";
auto dummy_dm = create_dummy_dm();
auto apply_restrict = [&]() {
auto get_is = [](
auto v) {
IS iy;
auto create = [&]() {
CHKERR VecGetOwnershipRange(
v, &ystart, NULL);
CHKERR ISCreateStride(PETSC_COMM_SELF,
n, ystart, 1, &iy);
};
};
auto iy = get_is(glob_x);
DMSubDomainRestrict(
simple->getDM(), s, PETSC_NULL, dummy_dm),
"restrict");
"get X0");
DMGetNamedGlobalVector(dummy_dm, "TSTheta_Xdot", &Xdot),
"get Xdot");
auto forward_ghost = [](
auto g) {
CHKERR VecGhostUpdateBegin(
g, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
g, INSERT_VALUES, SCATTER_FORWARD);
};
<< "Reverse restrict: X0 " << get_norm(X0) << " Xdot "
<< get_norm(Xdot);
}
return std::vector<Vec>{X0, Xdot};
};
auto ts_solver_vecs = apply_restrict();
if (ts_solver_vecs.size()) {
int ii = 0;
for (
auto v : ts_solver_vecs) {
MOFEM_LOG(
"FS", Sev::inform) <<
"Solve projection vector";
SCATTER_REVERSE);
for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
MOFEM_LOG(
"WORLD", Sev::verbose) <<
"Zero field " <<
f;
CHKERR field_blas->fieldLambdaOnEntities(zero_dofs,
f);
}
SCATTER_REVERSE);
SCATTER_FORWARD);
MOFEM_LOG(
"FS", Sev::inform) <<
"Norm V " << get_norm(
v);
}
CHKERR DMRestoreNamedGlobalVector(dummy_dm,
"TSTheta_X0",
&ts_solver_vecs[0]);
CHKERR DMRestoreNamedGlobalVector(dummy_dm,
"TSTheta_Xdot",
&ts_solver_vecs[1]);
}
for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
MOFEM_LOG(
"WORLD", Sev::verbose) <<
"Zero field " <<
f;
CHKERR field_blas->fieldLambdaOnEntities(zero_dofs,
f);
}
}
};
auto post_proc = [&](auto exe_test) {
auto post_proc_fe = boost::make_shared<PostProcEleDomain>(mField);
auto &pip = post_proc_fe->getOpPtrVector();
post_proc_fe->exeTestHook = exe_test;
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
fe_parent->getOpPtrVector(), {H1});
return fe_parent;
},
auto u_ptr = boost::make_shared<MatrixDouble>();
auto p_ptr = boost::make_shared<VectorDouble>();
auto h_ptr = boost::make_shared<VectorDouble>();
auto g_ptr = boost::make_shared<VectorDouble>();
CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"U",
CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"P",
CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"H",
CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"G",
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
{{"P", p_ptr}, {"H", h_ptr}, {"G", g_ptr}},
{{"U", u_ptr}},
{},
{}
)
);
CHKERR post_proc_fe->writeFile(
"out_projection.h5m");
};
});
});
}
fe_domain_lhs.swap(pip_mng->getDomainLhsFE());
fe_domain_rhs.swap(pip_mng->getDomainRhsFE());
fe_bdy_lhs.swap(pip_mng->getBoundaryLhsFE());
fe_bdy_rhs.swap(pip_mng->getBoundaryRhsFE());
}
auto add_parent_field_domain = [&](auto fe, auto op, auto field) {
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
};
auto add_parent_field_bdy = [&](auto fe, auto op, auto field) {
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
};
};
auto dot_u_ptr = boost::make_shared<MatrixDouble>();
auto u_ptr = boost::make_shared<MatrixDouble>();
auto grad_u_ptr = boost::make_shared<MatrixDouble>();
auto dot_h_ptr = boost::make_shared<VectorDouble>();
auto h_ptr = boost::make_shared<VectorDouble>();
auto grad_h_ptr = boost::make_shared<MatrixDouble>();
auto g_ptr = boost::make_shared<VectorDouble>();
auto grad_g_ptr = boost::make_shared<MatrixDouble>();
auto lambda_ptr = boost::make_shared<VectorDouble>();
auto p_ptr = boost::make_shared<VectorDouble>();
auto div_u_ptr = boost::make_shared<VectorDouble>();
auto set_domain_general = [&](auto fe) {
auto &pip = fe->getOpPtrVector();
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
fe_parent->getOpPtrVector(), {H1});
return fe_parent;
},
CHKERR add_parent_field_domain(fe, UDO::OPROW,
"U");
pip.push_back(
"U", grad_u_ptr));
pip.push_back(
"U", div_u_ptr));
break;
pip.push_back(
"U", div_u_ptr));
break;
default:
"Coordinate system not implemented");
}
CHKERR add_parent_field_domain(fe, UDO::OPROW,
"H");
pip.push_back(
CHKERR add_parent_field_domain(fe, UDO::OPROW,
"G");
pip.push_back(
CHKERR add_parent_field_domain(fe, UDO::OPROW,
"P");
};
auto set_domain_rhs = [&](auto fe) {
auto &pip = fe->getOpPtrVector();
CHKERR set_domain_general(fe);
CHKERR add_parent_field_domain(fe, UDO::OPROW,
"U");
pip.push_back(
new OpRhsU(
"U", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr,
grad_h_ptr, g_ptr, p_ptr));
CHKERR add_parent_field_domain(fe, UDO::OPROW,
"H");
pip.push_back(
new OpRhsH<false>(
"H", u_ptr, dot_h_ptr, h_ptr, grad_h_ptr,
grad_g_ptr));
CHKERR add_parent_field_domain(fe, UDO::OPROW,
"G");
CHKERR add_parent_field_domain(fe, UDO::OPROW,
"P");
"P", div_u_ptr, [](
const double r,
const double,
const double) {
}));
"P", p_ptr, [](
const double r,
const double,
const double) {
}));
};
auto set_domain_lhs = [&](auto fe) {
auto &pip = fe->getOpPtrVector();
CHKERR set_domain_general(fe);
CHKERR add_parent_field_domain(fe, UDO::OPROW,
"U");
{
CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"U");
pip.push_back(
new OpLhsU_dU(
"U", u_ptr, grad_u_ptr, h_ptr));
CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"H");
pip.push_back(
new OpLhsU_dH(
"U",
"H", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr, g_ptr));
CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"G");
pip.push_back(
new OpLhsU_dG(
"U",
"G", grad_h_ptr));
}
CHKERR add_parent_field_domain(fe, UDO::OPROW,
"H");
{
CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"U");
pip.push_back(
new OpLhsH_dU(
"H",
"U", grad_h_ptr));
CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"H");
CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"G");
}
CHKERR add_parent_field_domain(fe, UDO::OPROW,
"G");
{
CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"H");
CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"G");
}
CHKERR add_parent_field_domain(fe, UDO::OPROW,
"P");
{
CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"U");
"P", "U",
[](
const double r,
const double,
const double) {
},
true, false));
break;
"P", "U",
[](
const double r,
const double,
const double) {
},
true, false));
break;
default:
"Coordinate system not implemented");
}
CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"P");
pip.push_back(
new OpDomainMassP(
"P",
"P", [](
double r,
double,
double) {
}));
}
};
auto get_block_name = [](auto name) {
return boost::format("%s(.*)") % "WETTING_ANGLE";
};
auto get_blocks = [&](auto &&name) {
std::regex(name.str()));
};
auto set_boundary_rhs = [&](auto fe) {
auto &pip = fe->getOpPtrVector();
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"U");
CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"L");
pip,
mField,
"L", {},
"INFLUX",
CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"U");
auto wetting_block = get_blocks(get_block_name("WETTING_ANGLE"));
if (wetting_block.size()) {
auto op_bdy_side =
op_bdy_side->getOpPtrVector(), {H1});
op_bdy_side->getSideFEPtr(), "", UDO::OPSPACE,
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
fe_parent->getOpPtrVector(), {H1});
return fe_parent;
},
CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPROW,
"H");
op_bdy_side->getOpPtrVector().push_back(
pip.push_back(op_bdy_side);
for (auto &b : wetting_block) {
std::vector<double> attr_vec;
CHKERR b->getMeshsetIdEntitiesByDimension(
b->getAttributes(attr_vec);
if (attr_vec.size() != 1)
"Should be one attribute");
MOFEM_LOG(
"FS", Sev::inform) <<
"Wetting angle: " << attr_vec.front();
CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"G");
"G", grad_h_ptr, boost::make_shared<Range>(force_edges),
attr_vec.front()));
}
}
};
auto set_boundary_lhs = [&](auto fe) {
auto &pip = fe->getOpPtrVector();
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"L");
CHKERR add_parent_field_bdy(fe, UDO::OPCOL,
"U");
auto wetting_block = get_blocks(get_block_name("WETTING_ANGLE"));
if (wetting_block.size()) {
auto col_ind_ptr = boost::make_shared<std::vector<VectorInt>>();
auto col_diff_base_ptr = boost::make_shared<std::vector<MatrixDouble>>();
auto op_bdy_side =
op_bdy_side->getOpPtrVector(), {H1});
op_bdy_side->getSideFEPtr(), "", UDO::OPSPACE,
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
fe_parent->getOpPtrVector(), {H1});
return fe_parent;
},
CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPROW,
"H");
op_bdy_side->getOpPtrVector().push_back(
CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPCOL,
"H");
op_bdy_side->getOpPtrVector().push_back(
pip.push_back(op_bdy_side);
for (auto &b : wetting_block) {
std::vector<double> attr_vec;
CHKERR b->getMeshsetIdEntitiesByDimension(
b->getAttributes(attr_vec);
if (attr_vec.size() != 1)
"Should be one attribute");
<< "wetting angle edges size " << force_edges.size();
CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"G");
"G", grad_h_ptr, col_ind_ptr, col_diff_base_ptr,
boost::make_shared<Range>(force_edges), attr_vec.front()));
}
}
};
CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
CHKERR set_boundary_rhs(pip_mng->getCastBoundaryRhsFE());
CHKERR set_boundary_lhs(pip_mng->getCastBoundaryLhsFE());
}
boost::shared_ptr<PostProcEleDomainCont> post_proc,
boost::shared_ptr<PostProcEleBdyCont> post_proc_edge,
std::pair<boost::shared_ptr<BoundaryEle>, boost::shared_ptr<VectorDouble>>
p)
MOFEM_LOG(
"FS", Sev::verbose) <<
"Mesh pre proc";
auto post_proc_begin =
boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
CHKERR post_proc_end->writeFile(
"out_step_" + boost::lexical_cast<std::string>(
ts_step) +
".h5m");
MOFEM_LOG(
"FS", Sev::verbose) <<
"Mesh pre proc done";
}
MPI_COMM_WORLD);
<< " lift vec x: " << (*liftVec)[0] << " y: " << (*liftVec)[1];
}
private:
boost::shared_ptr<PostProcEleDomainCont>
postProc;
boost::shared_ptr<BoundaryEle>
liftFE;
boost::shared_ptr<VectorDouble>
liftVec;
};
DM subdm;
auto setup_subdm = [&](auto dm) {
auto level_ents_ptr = boost::make_shared<Range>();
CHKERR bit_mng->getEntitiesByRefLevel(
CHKERR DMSetType(subdm,
"DMMOFEM");
for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
}
};
};
auto get_fe_post_proc = [&](auto post_proc_mesh) {
auto add_parent_field_domain = [&](auto fe, auto op, auto field) {
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
};
auto post_proc_fe =
boost::make_shared<PostProcEleDomainCont>(
mField, post_proc_mesh);
post_proc_fe->exeTestHook = [](
FEMethod *fe_ptr) {
};
auto u_ptr = boost::make_shared<MatrixDouble>();
auto grad_u_ptr = boost::make_shared<MatrixDouble>();
auto h_ptr = boost::make_shared<VectorDouble>();
auto grad_h_ptr = boost::make_shared<MatrixDouble>();
auto p_ptr = boost::make_shared<VectorDouble>();
auto g_ptr = boost::make_shared<VectorDouble>();
auto grad_g_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector(), {H1});
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
fe_parent->getOpPtrVector(), {H1});
return fe_parent;
},
CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"U");
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
grad_u_ptr));
CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"H");
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"P");
post_proc_fe->getOpPtrVector().push_back(
CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"G");
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"H", h_ptr}, {"P", p_ptr}, {"G", g_ptr}},
{{"U", u_ptr}, {"H_GRAD", grad_h_ptr}, {"G_GRAD", grad_g_ptr}},
{{"GRAD_U", grad_u_ptr}},
{}
)
);
return post_proc_fe;
};
auto get_bdy_post_proc_fe = [&](auto post_proc_mesh) {
auto add_parent_field_bdy = [&](auto fe, auto op, auto field) {
return setParentDofs(
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
};
auto post_proc_fe =
boost::make_shared<PostProcEleBdyCont>(mField, post_proc_mesh);
post_proc_fe->exeTestHook = [](
FEMethod *fe_ptr) {
};
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
OpGetNormal(boost::shared_ptr<VectorDouble> l_ptr,
boost::shared_ptr<MatrixDouble> n_ptr)
ptrNormal(n_ptr) {}
auto t_n_fe = getFTensor1NormalsAtGaussPts();
ptrNormal->resize(
SPACE_DIM, getGaussPts().size2(),
false);
auto t_n = getFTensor1FromMat<SPACE_DIM>(*ptrNormal);
for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
t_n(
i) = t_n_fe(
i) * t_l / std::sqrt(t_n_fe(
i) * t_n_fe(
i));
++t_n_fe;
++t_l;
++t_n;
}
};
protected:
boost::shared_ptr<VectorDouble> ptrL;
boost::shared_ptr<MatrixDouble> ptrNormal;
};
auto u_ptr = boost::make_shared<MatrixDouble>();
auto p_ptr = boost::make_shared<VectorDouble>();
auto lambda_ptr = boost::make_shared<VectorDouble>();
auto normal_l_ptr = boost::make_shared<MatrixDouble>();
CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW,
"U");
post_proc_fe->getOpPtrVector().push_back(
CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW,
"L");
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpGetNormal(lambda_ptr, normal_l_ptr));
CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW,
"P");
post_proc_fe->getOpPtrVector().push_back(
op_ptr->doWorkRhsHook = [&](
DataOperator *base_op_ptr,
int side,
};
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
)
);
return post_proc_fe;
};
auto get_lift_fe = [&]() {
auto add_parent_field_bdy = [&](auto fe, auto op, auto field) {
return setParentDofs(
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
};
auto fe = boost::make_shared<BoundaryEle>(mField);
fe->exeTestHook = [](
FEMethod *fe_ptr) {
};
auto lift_ptr = boost::make_shared<VectorDouble>();
auto p_ptr = boost::make_shared<VectorDouble>();
auto ents_ptr = boost::make_shared<Range>();
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
std::vector<const CubitMeshSets *> vec_ptr;
std::regex("LIFT"), vec_ptr);
for (auto m_ptr : vec_ptr) {
ents, true);
ents_ptr->merge(ents);
}
MOFEM_LOG(
"FS", Sev::noisy) <<
"Lift ents " << (*ents_ptr);
CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"P");
fe->getOpPtrVector().push_back(
fe->getOpPtrVector().push_back(
return std::make_pair(fe, lift_ptr);
};
auto set_post_proc_monitor = [&](auto ts) {
DM dm;
boost::shared_ptr<FEMethod> null_fe;
auto post_proc_mesh = boost::make_shared<moab::Core>();
auto monitor_ptr = boost::make_shared<Monitor>(
get_fe_post_proc(post_proc_mesh), get_bdy_post_proc_fe(post_proc_mesh),
get_lift_fe());
null_fe, monitor_ptr);
};
auto ts_pre_post_proc = boost::make_shared<TSPrePostProc>();
ptr->fsRawPtr = this;
ptr->solverSubDM = create_solver_dm(
simple->getDM());
SCATTER_FORWARD);
CHKERR VecAssemblyBegin(ptr->globSol);
CHKERR VecAssemblyEnd(ptr->globSol);
auto sub_ts = pip_mng->createTSIM(ptr->solverSubDM);
CHKERR set_post_proc_monitor(sub_ts);
double ftime = 1;
auto print_fields_in_section = [&]() {
PetscInt num_fields;
CHKERR PetscSectionGetNumFields(section, &num_fields);
for (
int f = 0;
f < num_fields; ++
f) {
}
};
CHKERR print_fields_in_section();
CHKERR TSSolve(ts, ptr->globSol);
}
}
int main(
int argc,
char *argv[]) {
#ifdef PYTHON_INIT_SURFACE
Py_Initialize();
#endif
const char param_file[] = "param_file.petsc";
auto core_log = logging::core::get();
core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "FS"));
LogManager::setLog("FS");
try {
DMType dm_name = "DMMOFEM";
}
#ifdef PYTHON_INIT_SURFACE
if (Py_FinalizeEx() < 0) {
exit(120);
}
#endif
}
std::vector<Range>
"can not get vertices on bit 0");
Range plus_range, minus_range;
std::vector<EntityHandle> plus, minus;
for (auto p = vertices.pair_begin(); p != vertices.pair_end(); ++p) {
const auto s = p->second;
const auto lo_uid = DofEntity::getLoFieldEntityUId(field_bit_number,
f);
const auto hi_uid = DofEntity::getHiFieldEntityUId(field_bit_number, s);
auto it = dofs_mi.lower_bound(lo_uid);
const auto hi_it = dofs_mi.upper_bound(hi_uid);
plus.clear();
minus.clear();
plus.reserve(std::distance(it, hi_it));
minus.reserve(std::distance(it, hi_it));
for (; it != hi_it; ++it) {
const auto v = (*it)->getFieldData();
plus.push_back((*it)->getEnt());
else
minus.push_back((*it)->getEnt());
}
plus_range.insert_list(plus.begin(), plus.end());
minus_range.insert_list(minus.begin(), minus.end());
}
<< "Plus range " << plus_range << endl;
<< "Minus range " << minus_range << endl;
auto get_elems = [&](
auto &ents,
auto bit,
auto mask) {
moab::Interface::UNION),
"can not get adjacencies");
"can not filter elements with bit 0");
return adj;
};
CHKERR comm_mng->synchroniseEntities(plus_range);
CHKERR comm_mng->synchroniseEntities(minus_range);
auto common = intersect(ele_plus[0], ele_minus[0]);
ele_plus[0] = subtract(ele_plus[0], common);
ele_minus[0] = subtract(ele_minus[0], common);
auto get_children = [&](
auto &p,
auto &
c) {
CHKERR bit_mng->updateRangeByChildren(p,
c);
};
"get children");
"get children");
}
auto get_level = [&](
auto &p,
auto &
m,
auto z,
auto bit,
auto mask) {
"can not get vertices on bit");
for (
auto f = 0;
f != z; ++
f) {
l = get_elems(conn,
bit, mask);
}
};
vec_levels[
l] = get_level(ele_plus[
l], ele_minus[
l], 2 * overlap,
bit(
l),
}
std::string name = (boost::format(
"out_r%d.h5m") %
l).str();
"save mesh");
}
}
}
return vec_levels;
}
start_mask[s] = true;
CHKERR bit_mng->lambdaBitRefLevel(
true);
auto set_levels = [&](auto &&
vec_levels ) {
auto get_adj = [&](auto ents) {
"get conn");
ents.subset_by_dimension(
SPACE_DIM),
d,
false, conn,
moab::Interface::UNION),
"get adj");
}
ents.merge(conn);
return ents;
};
CHKERR bit_mng->updateRangeByParent(vec_levels[
l], parents);
level_prev = subtract(level_prev, parents);
level_prev.merge(vec_levels[
l]);
}
CHKERR bit_mng->getEntitiesByDimAndRefLevel(
level = get_adj(level);
}
};
auto set_skins = [&]() {
ParallelComm *pcomm =
"can't get bit level");
"can't get skin");
return bit_skin;
};
auto get_level_skin = [&]() {
CHKERR pcomm->filter_pstatus(skin_level_mesh,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, nullptr);
auto skin_level =
skin_level = subtract(skin_level,
skin_level_mesh);
true);
skin_level.merge(skin_level_verts);
level_prev);
skin.merge(subtract(skin_level, level_prev));
}
return skin;
};
auto resolve_shared = [&](auto &&skin) {
map<int, Range> map_procs;
tmp_skin, &map_procs);
for (
auto &
m : map_procs) {
from_other_procs.merge(
m.second);
}
}
auto common = intersect(
skin, from_other_procs);
skin.merge(from_other_procs);
if (!common.empty()) {
skin = subtract(skin, common);
}
return skin;
};
auto get_parent_level_skin = [&](auto skin) {
CHKERR bit_mng->updateRangeByParent(
skin.subset_by_dimension(
SPACE_DIM - 1), skin_parents);
true);
skin_parents.merge(skin_parent_verts);
skin_parents);
return skin_parents;
};
auto child_skin = resolve_shared(get_level_skin());
auto parent_skin = get_parent_level_skin(child_skin);
child_skin = subtract(child_skin, parent_skin);
};
auto set_current = [&]() {
last_level = subtract(last_level, skin_child);
};
std::string name = (boost::format(
"out_level%d.h5m") %
l).str();
"PARALLEL=WRITE_PART");
}
"MOAB", "PARALLEL=WRITE_PART");
"MOAB", "PARALLEL=WRITE_PART");
"MOAB", "PARALLEL=WRITE_PART");
"MOAB", "PARALLEL=WRITE_PART");
}
}
};
boost::shared_ptr<FEMethod> fe_top, std::string
field_name,
boost::function<boost::shared_ptr<ForcesAndSourcesCore>()> get_elem,
boost::function<void(boost::shared_ptr<ForcesAndSourcesCore>, int)>
add_parent_level =
[&](boost::shared_ptr<ForcesAndSourcesCore> parent_fe_pt, int level) {
if (level > 0) {
auto fe_ptr_current = get_elem();
add_parent_level(
boost::dynamic_pointer_cast<ForcesAndSourcesCore>(
fe_ptr_current),
level - 1);
if (op == ForcesAndSourcesCore::UserDataOperator::OPSPACE) {
parent_fe_pt->getOpPtrVector().push_back(
verbosity, sev));
} else {
parent_fe_pt->getOpPtrVector().push_back(
verbosity, sev));
}
}
};
add_parent_level(boost::dynamic_pointer_cast<ForcesAndSourcesCore>(fe_top),
}
auto cut_off_dofs = [&]() {
auto &m_field = ptr->fsRawPtr->mField;
CHKERR bit_mng->getEntitiesByTypeAndRefLevel(
auto cut_off_verts = [&](boost::shared_ptr<FieldEntity> ent_ptr) {
for (
auto &
h : ent_ptr->getEntFieldData()) {
}
};
CHKERR field_blas->fieldLambdaOnEntities(cut_off_verts,
"H",
¤t_verts);
};
}
MOFEM_LOG(
"FS", Sev::inform) <<
"Run step pre proc";
auto &m_field = ptr->fsRawPtr->mField;
auto get_norm = [&](auto x) {
double nrm;
CHKERR VecNorm(x, NORM_2, &nrm);
return nrm;
};
auto refine_problem = [&]() {
MOFEM_LOG(
"FS", Sev::inform) <<
"Refine problem";
CHKERR ptr->fsRawPtr->projectData();
};
auto set_jacobian_operators = [&]() {
};
auto set_solution = [&]() {
MOFEM_LOG(
"FS", Sev::inform) <<
"Set solution";
PetscObjectState state;
INSERT_VALUES, SCATTER_FORWARD);
<< "Set solution, vector norm " << get_norm(ptr->globSol);
};
PetscBool is_theta;
PetscObjectTypeCompare((PetscObject)ts, TSTHETA, &is_theta);
if (is_theta) {
CHKERR set_jacobian_operators();
} else {
"Sorry, only TSTheta handling is implemented");
}
PetscBarrier((PetscObject)ts);
}
}
auto &m_field = ptr->fsRawPtr->mField;
}
return 0;
}
auto sub_u = ptr->getSubVector();
auto scatter = ptr->getScatter(sub_u, u,
R);
auto apply_scatter_and_update = [&](auto x, auto sub_x) {
CHKERR VecScatterBegin(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
CHKERR VecScatterEnd(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateBegin(sub_x, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(sub_x, INSERT_VALUES, SCATTER_FORWARD);
};
CHKERR apply_scatter_and_update(u, sub_u);
CHKERR apply_scatter_and_update(u_t, sub_u_t);
CHKERR VecScatterBegin(scatter, sub_f,
f, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecScatterEnd(scatter, sub_f,
f, INSERT_VALUES, SCATTER_FORWARD);
}
}
PetscReal
a, Mat
A, Mat B,
void *ctx) {
auto sub_u = ptr->getSubVector();
auto scatter = ptr->getScatter(sub_u, u,
R);
auto apply_scatter_and_update = [&](auto x, auto sub_x) {
CHKERR VecScatterBegin(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
CHKERR VecScatterEnd(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateBegin(sub_x, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(sub_x, INSERT_VALUES, SCATTER_FORWARD);
};
CHKERR apply_scatter_and_update(u, sub_u);
CHKERR apply_scatter_and_update(u_t, sub_u_t);
ptr->tsCtxPtr.get());
}
}
auto get_norm = [&](auto x) {
double nrm;
CHKERR VecNorm(x, NORM_2, &nrm);
return nrm;
};
auto sub_u = ptr->getSubVector();
auto scatter = ptr->getScatter(sub_u, u,
R);
CHKERR VecScatterBegin(scatter, u, sub_u, INSERT_VALUES, SCATTER_REVERSE);
CHKERR VecScatterEnd(scatter, u, sub_u, INSERT_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateBegin(sub_u, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(sub_u, INSERT_VALUES, SCATTER_FORWARD);
<< "u norm " << get_norm(u) << " u sub nom " << get_norm(sub_u);
}
}
MOFEM_LOG(
"FS", Sev::verbose) <<
"SetUP sub PC";
ptr->subKSP =
createKSP(ptr->fsRawPtr->mField.get_comm());
CHKERR KSPSetFromOptions(ptr->subKSP);
}
};
auto sub_x = ptr->getSubVector();
auto scatter = ptr->getScatter(sub_x, pc_x,
R);
CHKERR VecScatterBegin(scatter, pc_f, sub_f, INSERT_VALUES,
SCATTER_REVERSE);
CHKERR VecScatterEnd(scatter, pc_f, sub_f, INSERT_VALUES, SCATTER_REVERSE);
CHKERR KSPSetOperators(ptr->subKSP, ptr->subB, ptr->subB);
MOFEM_LOG(
"FS", Sev::verbose) <<
"PCShell solve";
CHKERR KSPSolve(ptr->subKSP, sub_f, sub_x);
MOFEM_LOG(
"FS", Sev::verbose) <<
"PCShell solve <- done";
CHKERR VecScatterBegin(scatter, sub_x, pc_x, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR VecScatterEnd(scatter, sub_x, pc_x, INSERT_VALUES, SCATTER_FORWARD);
}
};
auto prb_ptr = ptr->fsRawPtr->mField.get_problem("SUB_SOLVER");
if (auto sub_data = prb_ptr->getSubData()) {
auto is = sub_data->getSmartColIs();
VecScatter s;
"crate scatter");
} else {
"crate scatter");
}
}
}
}
}
SNES snes;
auto set_section_monitor = [&](auto snes) {
(void *)(snes_ctx_ptr.get()), nullptr);
};
CHKERR set_section_monitor(snes);
CHKERR KSPSetType(ksp, KSPPREONLY);
CHKERR PCSetType(sub_pc, PCSHELL);
}