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
IntegrationType::GAUSS;
template <int DIM>
using SideOp = SideEle::UserDataOperator;
template <CoordinateTypes COORD_TYPE>
BoundaryNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 1>;
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;
};
auto save_range = [](moab::Interface &moab,
const std::string name,
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;
auto lo_it = dofs->lower_bound(FieldEntity::getLoBitNumberUId(bit_number));
auto hi_it = dofs->upper_bound(FieldEntity::getHiBitNumberUId(bit_number));
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());
for (auto d : *dofs) {
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:
Vec f,
void *ctx);
PetscReal a, Mat
A, Mat
B,
void *ctx);
void *ctx);
boost::shared_ptr<SnesCtx>
boost::shared_ptr<TsCtx>
};
private:
boost::shared_ptr<FEMethod> fe_top, std::string
field_name,
ForcesAndSourcesCore::UserDataOperator::OpType op,
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) {
[&]() {
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 *))MoFEMSNESMonitorFields,
(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) {
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
};
auto add_parent_field_bdy = [&](
auto fe,
auto op,
auto field,
auto bit) {
[&]() {
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;
},
};
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
};
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");
CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpFluidFlux>::add(
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>>
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) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
};
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) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
};
[&]() {
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 = getFTensor1FromMat<SPACE_DIM>(*ptrNormal);
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(),
OpPPMap::DataMapVec{{"P", p_ptr}},
OpPPMap::DataMapMat{{"U", u_ptr}, {"L", normal_l_ptr}},
OpPPMap::DataMapMat(),
OpPPMap::DataMapMat()
)
);
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) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
};
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
auto core_log = logging::core::get();
core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "FS"));
LogManager::setLog("FS");
try {
DMType dm_name = "DMMOFEM";
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
}
#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,
ForcesAndSourcesCore::UserDataOperator::OpType op,
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;
}
Vec f, 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);
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());
}
}
Vec u, void *ctx) {
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 *))MoFEMSNESMonitorFields,
(void *)(snes_ctx_ptr.get()), nullptr);
};
CHKERR set_section_monitor(snes);
CHKERR KSPSetType(ksp, KSPPREONLY);
CHKERR PCSetType(sub_pc, PCSHELL);
}
#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
static boost::weak_ptr< TSPrePostProc > tsPrePostProc
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
double lambda
surface tension
PostProcBrokenMeshInMoabBaseCont< DomainEle > PostProcEleDomainCont
OpDomainMassH OpDomainMassG
auto init_h
Initialisation function.
FTensor::Index< 'l', SPACE_DIM > l
PostProcBrokenMeshInMoabBaseCont< BoundaryEle > PostProcEleBdyCont
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
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
BoundaryNaturalBC::OpFlux< NaturalMeshsetType< BLOCKSET >, 1, 1 > OpFluidFlux
auto get_current_bit
dofs bit used to do calculations
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
#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 marker
set bit to marker
auto test_bit_child
lambda function used to select elements on which finite element pipelines are executed.
FTensor::Index< 'i', SPACE_DIM > i
constexpr CoordinateTypes coord_type
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
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const FTensor::Tensor2< T, Dim, Dim > Vec
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
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)
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.
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
constexpr IntegrationType I
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]
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.
Data on single entity (This is passed as argument to DataOperator::doWork)
default operator for TRI element
auto getFTensor1NormalsAtGaussPts()
get normal at integration points
@ 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.
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 rate of scalar field at integration points.
Get value at integration points for scalar 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.
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
boost::shared_ptr< VectorDouble > liftVec
boost::shared_ptr< BoundaryEle > liftFE
boost::shared_ptr< PostProcEle > postProc
boost::shared_ptr< PostProcEleBdyCont > postProcEdge
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
NaturalBC< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS > BoundaryNaturalBC
[Body and heat source]