v0.14.0
Loading...
Searching...
No Matches
free_surface.cpp

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.

Implementation based on [49]

/**
* \file free_surface.cpp
* \example free_surface.cpp
*
* 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.
*
* Implementation based on \cite Lovric2019-qn
*/
#include <MoFEM.hpp>
#include <petsc/private/petscimpl.h>
using namespace MoFEM;
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;
MoFEMErrorCode surfaceInit(const std::string py_file) {
try {
// create main module
auto main_module = bp::import("__main__");
mainNamespace = main_module.attr("__dict__");
bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
// create a reference to python function
surfaceFun = mainNamespace["surface"];
} catch (bp::error_already_set const &) {
// print all other errors to stderr
PyErr_Print();
}
};
MoFEMErrorCode evalSurface(
double x, double y, double z, double eta, double &s
) {
try {
// call python function
s = bp::extract<double>(surfaceFun(x, y, z, eta));
} catch (bp::error_already_set const &) {
// print all other errors to stderr
PyErr_Print();
}
}
private:
bp::object mainNamespace;
bp::object surfaceFun;
};
static boost::weak_ptr<SurfacePython> surfacePythonWeakPtr;
#endif
int coord_type = EXECUTABLE_COORD_TYPE;
constexpr int BASE_DIM = 1;
constexpr int SPACE_DIM = 2;
constexpr int U_FIELD_DIM = SPACE_DIM;
constexpr AssemblyType A = AssemblyType::PETSC; //< selected assembly type
constexpr IntegrationType I =
IntegrationType::GAUSS; //< selected integration type
template <int DIM>
using SideOp = SideEle::UserDataOperator;
template <CoordinateTypes COORD_TYPE>
// Flux is applied by Lagrange Multiplie
using OpFluidFlux =
BoundaryNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 1>;
// mesh refinement
int order = 3; ///< approximation order
int nb_levels = 4; //< number of refinement levels
int refine_overlap = 4; //< mesh overlap while refine
constexpr bool debug = true;
auto get_start_bit = []() {
return nb_levels + 1;
}; //< first refinement level for computational mesh
auto get_current_bit = []() {
return 2 * get_start_bit() + 1;
}; ///< dofs bit used to do calculations
auto get_skin_parent_bit = []() { return 2 * get_start_bit() + 2; };
auto get_skin_child_bit = []() { return 2 * get_start_bit() + 3; };
auto get_projection_bit = []() { return 2 * get_start_bit() + 4; };
auto get_skin_projection_bit = []() { return 2 * get_start_bit() + 5; };
// FIXME: Set parameters from command line
// Physical parameters
double a0 = 980;
double rho_m = 0.998;
double mu_m = 0.010101 * 1e2;
double rho_p = 0.0012;
double mu_p = 0.000182 * 1e2;
double lambda = 73; ///< surface tension
double W = 0.25;
// Model parameters
double h = 0.01 / 8; // mesh size
double eta = h;
double eta2 = eta * eta;
// Numerical parameters
double md = 1e-2; // mobility
double eps = 1e-12;
double tol = std::numeric_limits<float>::epsilon();
double rho_ave = (rho_p + rho_m) / 2;
double rho_diff = (rho_p - rho_m) / 2;
double mu_ave = (mu_p + mu_m) / 2;
double mu_diff = (mu_p - mu_m) / 2;
double kappa = (3. / (4. * std::sqrt(2. * W))) * (lambda / eta);
auto integration_rule = [](int, int, int) { return 2 * order + 1; };
auto cylindrical = [](const double r) {
return 2 * M_PI * r;
else
return 1.;
};
auto wetting_angle_sub_stepping = [](auto ts_step) {
constexpr int sub_stepping = 16;
return std::min(1., static_cast<double>(ts_step) / sub_stepping);
};
// cut off function
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; };
auto cut_off = [](const double h) { return my_max(my_min(h)); };
auto d_cut_off = [](const double h) {
if (h >= -1 && h < 1)
return 1.;
else
return 0.;
};
auto phase_function = [](const double h, const double diff, const double ave) {
return diff * h + ave;
};
auto d_phase_function_h = [](const double h, const double diff) {
return diff;
};
// order function (free energy)
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); };
auto get_M0 = [](auto h) { return md; };
auto get_M0_dh = [](auto h) { return 0; };
auto get_M2 = [](auto h) {
return md * (1 - h * h);
};
auto get_M2_dh = [](auto h) { return -md * 2 * h; };
auto get_M3 = [](auto h) {
const double h2 = h * h;
const double h3 = h2 * h;
if (h >= 0)
return md * (2 * h3 - 3 * h2 + 1);
else
return md * (-2 * h3 - 3 * h2 + 1);
};
auto get_M3_dh = [](auto h) {
if (h >= 0)
return md * (6 * h * (h - 1));
else
return md * (-6 * h * (h + 1));
};
auto get_M = [](auto h) { return get_M0(h); };
auto get_M_dh = [](auto h) { return get_M0_dh(h); };
auto get_D = [](const double A) {
t_D(i, j, k, l) = A * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
return t_D;
};
// some examples of initialisation functions
auto kernel_oscillation = [](double r, double y, double) {
constexpr int n = 3;
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)));
};
auto kernel_eye = [](double r, double y, double) {
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)));
};
auto capillary_tube = [](double x, double y, double z) {
constexpr double water_height = 0.;
return tanh((water_height - y) / (eta * std::sqrt(2)));
;
};
auto bubble_device = [](double x, double y, double z) {
return -tanh((-0.039 - x) / (eta * std::sqrt(2)));
};
/**
* @brief Initialisation function
*
* @note If UMs are compiled with Python to initialise phase field "H"
* surface.py function is used, which has to be present in execution folder.
*
*/
auto init_h = [](double r, double y, double theta) {
#ifdef PYTHON_INIT_SURFACE
double s = 1;
if (auto ptr = surfacePythonWeakPtr.lock()) {
CHK_THROW_MESSAGE(ptr->evalSurface(r, y, theta, eta, s),
"error eval python");
}
return s;
#else
return bubble_device(r, y, theta);
// return capillary_tube(r, y, theta);
// return kernel_eye(r, y, theta);
#endif
};
auto wetting_angle = [](double water_level) { return water_level; };
auto bit = [](auto b) { return BitRefLevel().set(b); };
auto marker = [](auto b) { return BitRefLevel().set(BITREFLEVEL_SIZE - b); };
auto get_fe_bit = [](FEMethod *fe_ptr) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel();
};
auto get_global_size = [](int l_size) {
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,
const Range r) {
if (get_global_size(r.size())) {
auto out_meshset = get_temp_meshset_ptr(moab);
CHKERR moab.add_entities(*out_meshset, r);
CHKERR moab.write_file(name.c_str(), "MOAB", "PARALLEL=WRITE_PART",
out_meshset->get_ptr(), 1);
}
};
/**
* @brief get entities of dofs in the problem - used for debugging
*
*/
auto get_dofs_ents_by_field_name = [](auto dm, auto field_name) {
auto prb_ptr = getProblemPtr(dm);
std::vector<EntityHandle> ents_vec;
MoFEM::Interface *m_field_ptr;
CHKERR DMoFEMGetInterfacePtr(dm, &m_field_ptr);
auto bit_number = m_field_ptr->get_field_bit_number(field_name);
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);
return r;
};
/**
* @brief get entities of dofs in the problem - used for debugging
*
*/
auto get_dofs_ents_all = [](auto dm) {
auto prb_ptr = getProblemPtr(dm);
std::vector<EntityHandle> ents_vec;
MoFEM::Interface *m_field_ptr;
CHKERR DMoFEMGetInterfacePtr(dm, &m_field_ptr);
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);
return r;
};
using namespace FreeSurfaceOps;
struct FreeSurface;
enum FR { F, R }; // F - forward, and reverse
/**
* @brief Set of functions called by PETSc solver used to refine and update
* mesh.
*
* @note Currently theta method is only handled by this code.
*
*/
struct TSPrePostProc {
TSPrePostProc() = default;
virtual ~TSPrePostProc() = default;
/**
* @brief Used to setup TS solver
*
* @param ts
* @return MoFEMErrorCode
*/
SmartPetscObj<VecScatter> getScatter(Vec x, Vec y, enum FR fr);
private:
/**
* @brief Pre process time step
*
* Refine mesh and update fields
*
* @param ts
* @return MoFEMErrorCode
*/
static MoFEMErrorCode tsPreProc(TS ts);
/**
* @brief Post process time step.
*
* Currently that function do not make anything major
*
* @param ts
* @return MoFEMErrorCode
*/
static MoFEMErrorCode tsPostProc(TS ts);
static MoFEMErrorCode tsPreStage(TS ts);
static MoFEMErrorCode tsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t,
Vec f,
void *ctx); //< Wrapper for SNES Rhs
static MoFEMErrorCode tsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t,
PetscReal a, Mat A, Mat B,
void *ctx); ///< Wrapper for SNES Lhs
static MoFEMErrorCode tsMonitor(TS ts, PetscInt step, PetscReal t, Vec u,
void *ctx); ///< Wrapper for TS monitor
static MoFEMErrorCode pcSetup(PC pc);
static MoFEMErrorCode pcApply(PC pc, Vec pc_f, Vec pc_x);
SmartPetscObj<Vec> globRes; //< global residual
SmartPetscObj<Mat> subB; //< sub problem tangent matrix
SmartPetscObj<KSP> subKSP; //< sub problem KSP solver
boost::shared_ptr<SnesCtx>
snesCtxPtr; //< infernal data (context) for MoFEM SNES fuctions
boost::shared_ptr<TsCtx>
tsCtxPtr; //< internal data (context) for MoFEM TS functions.
};
static boost::weak_ptr<TSPrePostProc> tsPrePostProc;
struct FreeSurface {
FreeSurface(MoFEM::Interface &m_field) : mField(m_field) {}
private:
/// @brief Find entities on refinement levels
/// @param overlap level of overlap
/// @return
std::vector<Range> findEntitiesCrossedByPhaseInterface(size_t overlap);
/// @brief
/// @param ents
/// @param level
/// @param mask
/// @return
Range findParentsToRefine(Range ents, BitRefLevel level, BitRefLevel mask);
/// @brief
/// @param overlap
/// @return
MoFEMErrorCode refineMesh(size_t overlap);
/// @brief Create hierarchy of elements run on parents levels
/// @param fe_top pipeline element to which hierarchy is attached
/// @param op type of operator OPSPACE, OPROW, OPCOL or OPROWCOL
/// @param get_elema lambda function to create element on hierarchy
/// @param verbosity verbosity level
/// @param sev severity level
/// @return
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);
friend struct TSPrePostProc;
};
//! [Run programme]
}
//! [Run programme]
//! [Read mesh]
MOFEM_LOG("FS", Sev::inform) << "Read mesh for problem";
simple->getParentAdjacencies() = true;
simple->getBitRefLevel() = BitRefLevel();
CHKERR simple->getOptions();
CHKERR simple->loadFile();
}
//! [Read mesh]
//! [Set up problem]
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-nb_levels", &nb_levels,
PETSC_NULL);
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-refine_overlap", &refine_overlap,
PETSC_NULL);
const char *coord_type_names[] = {"cartesian", "polar", "cylindrical",
"spherical"};
CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-coords", coord_type_names,
MOFEM_LOG("FS", Sev::inform) << "Approximation order = " << order;
MOFEM_LOG("FS", Sev::inform)
<< "Number of refinement levels nb_levels = " << nb_levels;
nb_levels += 1;
auto bit_mng = mField.getInterface<BitRefManager>();
CHKERR PetscOptionsGetScalar(PETSC_NULL, "Acceleration", "-a0", &a0, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "\"Minus\" phase density rho_m",
"-rho_m", &rho_m, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "\"Minus\" phase viscosity", "-mu_m",
&mu_m, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "\"Plus\" phase density", "-rho_p",
&rho_p, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "\"Plus\" phase viscosity", "-mu_p",
&mu_p, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "Surface tension", "-lambda",
&lambda, PETSC_NULL);
"Height of the well in energy functional", "-W",
&W, PETSC_NULL);
rho_ave = (rho_p + rho_m) / 2;
rho_diff = (rho_p - rho_m) / 2;
mu_ave = (mu_p + mu_m) / 2;
mu_diff = (mu_p - mu_m) / 2;
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-h", &h, PETSC_NULL);
eta = h;
eta2 = eta * eta;
kappa = (3. / (4. * std::sqrt(2. * W))) * (lambda / eta);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-md", &eta, PETSC_NULL);
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;
MOFEM_LOG("FS", Sev::inform) << "Surface tension lambda = " << lambda;
MOFEM_LOG("FS", Sev::inform)
<< "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;
MOFEM_LOG("FS", Sev::inform) << "Average density rho_ave = " << rho_ave;
MOFEM_LOG("FS", Sev::inform) << "Difference density rho_diff = " << rho_diff;
MOFEM_LOG("FS", Sev::inform) << "Average viscosity mu_ave = " << mu_ave;
MOFEM_LOG("FS", Sev::inform) << "Difference viscosity mu_diff = " << mu_diff;
MOFEM_LOG("FS", Sev::inform) << "kappa = " << kappa;
// Fields on domain
// Velocity field
// Pressure field
CHKERR simple->addDomainField("P", H1, AINSWORTH_LEGENDRE_BASE, 1);
// Order/phase fild
CHKERR simple->addDomainField("H", H1, AINSWORTH_LEGENDRE_BASE, 1);
// Chemical potential
CHKERR simple->addDomainField("G", H1, AINSWORTH_LEGENDRE_BASE, 1);
// Field on boundary
CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE,
CHKERR simple->addBoundaryField("H", H1, AINSWORTH_LEGENDRE_BASE, 1);
CHKERR simple->addBoundaryField("G", H1, AINSWORTH_LEGENDRE_BASE, 1);
// Lagrange multiplier which constrains slip conditions
CHKERR simple->addBoundaryField("L", H1, AINSWORTH_LEGENDRE_BASE, 1);
CHKERR simple->setFieldOrder("U", order);
CHKERR simple->setFieldOrder("P", order - 1);
CHKERR simple->setFieldOrder("H", order);
CHKERR simple->setFieldOrder("G", order);
CHKERR simple->setFieldOrder("L", order);
// Initialise bit ref levels
auto set_problem_bit = [&]() {
// Set bits to build adjacencies between parents and children. That is
// used by simple interface.
simple->getBitAdjEnt() = BitRefLevel().set();
simple->getBitAdjParent() = BitRefLevel().set();
simple->getBitRefLevel() = BitRefLevel().set();
simple->getBitRefLevelMask() = BitRefLevel().set();
};
CHKERR set_problem_bit();
CHKERR simple->setUp();
}
//! [Set up problem]
//! [Boundary condition]
#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 pip_mng = mField.getInterface<PipelineManager>();
auto bc_mng = mField.getInterface<BcManager>();
auto bit_mng = mField.getInterface<BitRefManager>();
auto dm = simple->getDM();
using UDO = ForcesAndSourcesCore::UserDataOperator;
auto reset_bits = [&]() {
BitRefLevel start_mask;
for (auto s = 0; s != get_start_bit(); ++s)
start_mask[s] = true;
// reset bit ref levels
CHKERR bit_mng->lambdaBitRefLevel(
[&](EntityHandle ent, BitRefLevel &bit) { bit &= start_mask; });
Range level0;
CHKERR bit_mng->getEntitiesByRefLevel(bit(0), BitRefLevel().set(), level0);
CHKERR bit_mng->setNthBitRefLevel(level0, get_current_bit(), true);
CHKERR bit_mng->setNthBitRefLevel(level0, get_projection_bit(), true);
};
auto add_parent_field = [&](auto fe, auto op, auto field) {
return setParentDofs(
fe, field, op, bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
QUIET, Sev::noisy);
};
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();
fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
fe_parent->getOpPtrVector(), {H1});
return fe_parent;
},
QUIET, Sev::noisy);
CHKERR add_parent_field(fe, UDO::OPROW, "H");
pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
pip.push_back(
CHKERR add_parent_field(fe, UDO::OPROW, "G");
pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
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 DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
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) {
CHKERR set_generic(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) {
CHKERR set_generic(fe);
auto &pip = fe->getOpPtrVector();
CHKERR add_parent_field(fe, UDO::OPROW, "H");
CHKERR add_parent_field(fe, UDO::OPCOL, "H");
pip.push_back(new OpLhsH_dH<true>("H", nullptr, h_ptr, grad_g_ptr));
CHKERR add_parent_field(fe, UDO::OPCOL, "G");
pip.push_back(new OpLhsH_dG<true>("H", "G", h_ptr));
CHKERR add_parent_field(fe, UDO::OPROW, "G");
pip.push_back(new OpLhsG_dG("G"));
CHKERR add_parent_field(fe, UDO::OPCOL, "H");
pip.push_back(new OpLhsG_dH<true>("G", "H", h_ptr));
};
auto create_subdm = [&]() {
auto level_ents_ptr = boost::make_shared<Range>();
CHKERR mField.getInterface<BitRefManager>()->getEntitiesByRefLevel(
bit(get_current_bit()), BitRefLevel().set(), *level_ents_ptr);
DM subdm;
CHKERR DMCreate(mField.get_comm(), &subdm);
CHKERR DMSetType(subdm, "DMMOFEM");
CHKERR DMMoFEMCreateSubDM(subdm, dm, "SUB_INIT");
CHKERR DMMoFEMAddElement(subdm, simple->getDomainFEName());
CHKERR DMMoFEMSetSquareProblem(subdm, PETSC_TRUE);
CHKERR DMMoFEMSetDestroyProblem(subdm, PETSC_TRUE);
for (auto f : {"H", "G"}) {
CHKERR DMMoFEMAddSubFieldRow(subdm, f, level_ents_ptr);
CHKERR DMMoFEMAddSubFieldCol(subdm, f, level_ents_ptr);
}
CHKERR DMSetUp(subdm);
if constexpr (debug) {
if (mField.get_comm_size() == 1) {
auto dm_ents = get_dofs_ents_all(SmartPetscObj<DM>(subdm, true));
CHKERR save_range(mField.get_moab(), "sub_dm_init_ents_verts.h5m",
dm_ents.subset_by_type(MBVERTEX));
dm_ents = subtract(dm_ents, dm_ents.subset_by_type(MBVERTEX));
CHKERR save_range(mField.get_moab(), "sub_dm_init_ents.h5m", dm_ents);
}
}
return SmartPetscObj<DM>(subdm);
};
auto subdm = create_subdm();
pip_mng->getDomainRhsFE().reset();
pip_mng->getDomainLhsFE().reset();
CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule);
CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule);
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 D = createDMVector(subdm);
auto snes = pip_mng->createSNES(subdm);
auto snes_ctx_ptr = getDMSnesCtx(subdm);
auto set_section_monitor = [&](auto solver) {
CHKERR SNESMonitorSet(solver,
(MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
(void *)(snes_ctx_ptr.get()), nullptr);
};
auto print_section_field = [&]() {
auto section =
mField.getInterface<ISManager>()->sectionCreate("SUB_INIT");
PetscInt num_fields;
CHKERR PetscSectionGetNumFields(section, &num_fields);
for (int f = 0; f < num_fields; ++f) {
const char *field_name;
CHKERR PetscSectionGetFieldName(section, f, &field_name);
MOFEM_LOG("FS", Sev::inform)
<< "Field " << f << " " << std::string(field_name);
}
};
CHKERR set_section_monitor(snes);
CHKERR print_section_field();
for (auto f : {"U", "P", "H", "G", "L"}) {
CHKERR mField.getInterface<FieldBlas>()->setField(0, f);
}
CHKERR SNESSetFromOptions(snes);
CHKERR SNESSolve(snes, PETSC_NULL, D);
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(subdm, D, INSERT_VALUES, SCATTER_REVERSE);
};
CHKERR reset_bits();
CHKERR solve_init(
[](FEMethod *fe_ptr) { return get_fe_bit(fe_ptr).test(0); });
CHKERR refineMesh(refine_overlap);
for (auto f : {"U", "P", "H", "G", "L"}) {
CHKERR mField.getInterface<FieldBlas>()->setField(0, f);
}
CHKERR solve_init([](FEMethod *fe_ptr) {
return get_fe_bit(fe_ptr).test(get_start_bit() + nb_levels - 1);
});
CHKERR post_proc([](FEMethod *fe_ptr) {
return get_fe_bit(fe_ptr).test(get_start_bit() + nb_levels - 1);
});
pip_mng->getOpDomainRhsPipeline().clear();
pip_mng->getOpDomainLhsPipeline().clear();
// Remove DOFs where boundary conditions are set
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",
0, SPACE_DIM);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX", "L",
0, 0);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "ZERO",
"L", 0, 0);
}
//! [Boundary condition]
//! [Data projection]
// FIXME: Functionality in this method should be improved (projection should
// be done field by field), generalised, and move to become core
// functionality.
auto pip_mng = mField.getInterface<PipelineManager>();
auto bc_mng = mField.getInterface<BcManager>();
auto bit_mng = mField.getInterface<BitRefManager>();
auto field_blas = mField.getInterface<FieldBlas>();
// Store all existing elements pipelines, replace them by data projection
// pipelines, to put them back when projection is done.
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();
using UDO = ForcesAndSourcesCore::UserDataOperator;
// extract field data for domain parent element
auto add_parent_field_domain = [&](auto fe, auto op, auto field, auto bit) {
return setParentDofs(
fe, field, op, bit,
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
QUIET, Sev::noisy);
};
// extract field data for boundary parent element
auto add_parent_field_bdy = [&](auto fe, auto op, auto field, auto bit) {
return setParentDofs(
fe, field, op, bit,
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
QUIET, Sev::noisy);
};
// run this element on element with given bit, otherwise run other nested
// element
auto create_run_parent_op = [&](auto parent_fe_ptr, auto this_fe_ptr,
auto fe_bit) {
auto parent_mask = fe_bit;
parent_mask.flip();
return new OpRunParent(parent_fe_ptr, BitRefLevel().set(), parent_mask,
this_fe_ptr, fe_bit, BitRefLevel().set(), QUIET,
Sev::inform);
};
// create hierarchy of nested elements
auto get_parents_vel_fe_ptr = [&](auto this_fe_ptr, auto fe_bit) {
std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
for (int l = 0; l < nb_levels; ++l)
parents_elems_ptr_vec.emplace_back(
boost::make_shared<DomianParentEle>(mField));
for (auto l = 1; l < nb_levels; ++l) {
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];
};
// solve projection problem
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();
eval_fe_ptr, "", UDO::OPSPACE, bit(get_skin_projection_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
QUIET, Sev::noisy);
// That can be done much smarter, by block, field by field. For
// simplicity is like that.
CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW, "U",
pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW, "P",
pip.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW, "H",
pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW, "G",
pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
}
auto parent_eval_fe_ptr =
get_parents_vel_fe_ptr(eval_fe_ptr, bit(get_projection_bit()));
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();
assemble_fe_ptr, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
QUIET, Sev::noisy);
CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "U",
pip.push_back(new OpDomainAssembleVector("U", u_ptr));
CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "P",
pip.push_back(new OpDomainAssembleScalar("P", p_ptr));
CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "H",
pip.push_back(new OpDomainAssembleScalar("H", h_ptr));
CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "G",
pip.push_back(new OpDomainAssembleScalar("G", g_ptr));
}
auto parent_assemble_fe_ptr =
get_parents_vel_fe_ptr(assemble_fe_ptr, bit(get_current_bit()));
pip.push_back(create_run_parent_op(
parent_assemble_fe_ptr, assemble_fe_ptr, bit(get_current_bit())));
};
auto set_domain_lhs = [&](auto fe) {
auto &pip = fe->getOpPtrVector();
fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
QUIET, Sev::noisy);
// That can be done much smarter, by block, field by field. For simplicity
// is like that.
CHKERR add_parent_field_domain(fe, UDO::OPROW, "U",
CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U",
pip.push_back(new OpDomainMassU("U", "U"));
CHKERR add_parent_field_domain(fe, UDO::OPROW, "P",
CHKERR add_parent_field_domain(fe, UDO::OPCOL, "P",
pip.push_back(new OpDomainMassP("P", "P"));
CHKERR add_parent_field_domain(fe, UDO::OPROW, "H",
CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H",
pip.push_back(new OpDomainMassH("H", "H"));
CHKERR add_parent_field_domain(fe, UDO::OPROW, "G",
CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G",
pip.push_back(new OpDomainMassG("G", "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();
eval_fe_ptr, "", UDO::OPSPACE, bit(get_skin_projection_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
QUIET, Sev::noisy);
// That can be done much smarter, by block, field by field. For
// simplicity is like that.
CHKERR add_parent_field_bdy(eval_fe_ptr, UDO::OPROW, "L",
pip.push_back(new OpCalculateScalarFieldValues("L", l_ptr));
}
auto parent_eval_fe_ptr =
get_parents_vel_fe_ptr(eval_fe_ptr, bit(get_projection_bit()));
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();
assemble_fe_ptr, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
QUIET, Sev::noisy);
struct OpLSize : public BoundaryEleOp {
OpLSize(boost::shared_ptr<VectorDouble> l_ptr)
: BoundaryEleOp(NOSPACE, DomainEleOp::OPSPACE), lPtr(l_ptr) {}
MoFEMErrorCode doWork(int, EntityType, EntData &) {
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",
pip.push_back(new OpBoundaryAssembleScalar("L", l_ptr));
}
auto parent_assemble_fe_ptr =
get_parents_vel_fe_ptr(assemble_fe_ptr, bit(get_current_bit()));
pip.push_back(create_run_parent_op(
parent_assemble_fe_ptr, assemble_fe_ptr, bit(get_current_bit())));
};
auto set_bdy_lhs = [&](auto fe) {
auto &pip = fe->getOpPtrVector();
fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
QUIET, Sev::noisy);
// That can be done much smarter, by block, field by field. For simplicity
// is like that.
CHKERR add_parent_field_bdy(fe, UDO::OPROW, "L",
CHKERR add_parent_field_bdy(fe, UDO::OPCOL, "L",
pip.push_back(new OpBoundaryMassL("L", "L"));
};
auto create_subdm = [&]() -> SmartPetscObj<DM> {
auto level_ents_ptr = boost::make_shared<Range>();
CHKERR mField.getInterface<BitRefManager>()->getEntitiesByRefLevel(
bit(get_current_bit()), BitRefLevel().set(), *level_ents_ptr);
auto get_prj_ents = [&]() {
Range prj_mesh;
CHKERR bit_mng->getEntitiesByDimAndRefLevel(bit(get_projection_bit()),
BitRefLevel().set(),
SPACE_DIM, prj_mesh);
auto common_ents = intersect(prj_mesh, *level_ents_ptr);
prj_mesh = subtract(unite(*level_ents_ptr, prj_mesh), common_ents)
.subset_by_dimension(SPACE_DIM);
return prj_mesh;
};
auto prj_ents = get_prj_ents();
if (get_global_size(prj_ents.size())) {
auto rebuild = [&]() {
auto prb_mng = mField.getInterface<ProblemsManager>();
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);
// partition problem
CHKERR prb_mng->partitionFiniteElements("SUB_SOLVER", true, 0,
// set ghost nodes
CHKERR prb_mng->partitionGhostDofsOnDistributedMesh("SUB_SOLVER");
};
MOFEM_LOG("FS", Sev::verbose) << "Create projection problem";
CHKERR rebuild();
auto dm = simple->getDM();
DM subdm;
CHKERR DMCreate(mField.get_comm(), &subdm);
CHKERR DMSetType(subdm, "DMMOFEM");
CHKERR DMMoFEMCreateSubDM(subdm, dm, "SUB_SOLVER");
return SmartPetscObj<DM>(subdm);
}
MOFEM_LOG("FS", Sev::inform) << "Nothing to project";
};
auto create_dummy_dm = [&]() {
auto dummy_dm = createDM(mField.get_comm(), "DMMOFEM");
simple->getProblemName().c_str(),
BitRefLevel(), BitRefLevel()),
"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();
CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule);
CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule);
CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule);
CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule);
pip_mng->getDomainLhsFE()->exeTestHook = exe_test;
pip_mng->getDomainRhsFE()->exeTestHook = [](FEMethod *fe_ptr) {
return get_fe_bit(fe_ptr).test(nb_levels - 1);
};
pip_mng->getBoundaryLhsFE()->exeTestHook = exe_test;
pip_mng->getBoundaryRhsFE()->exeTestHook = [](FEMethod *fe_ptr) {
return get_fe_bit(fe_ptr).test(nb_levels - 1);
};
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 D = createDMVector(subdm);
auto F = vectorDuplicate(D);
auto ksp = pip_mng->createKSP(subdm);
CHKERR KSPSetFromOptions(ksp);
CHKERR KSPSetUp(ksp);
// get vector norm
auto get_norm = [&](auto x) {
double nrm;
CHKERR VecNorm(x, NORM_2, &nrm);
return nrm;
};
/**
* @brief Zero DOFs, used by FieldBlas
*/
auto zero_dofs = [](boost::shared_ptr<FieldEntity> ent_ptr) {
for (auto &v : ent_ptr->getEntFieldData()) {
v = 0;
}
};
auto solve = [&](auto S) {
CHKERR VecZeroEntries(S);
CHKERR VecZeroEntries(F);
CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
CHKERR KSPSolve(ksp, F, S);
CHKERR VecGhostUpdateBegin(S, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(S, INSERT_VALUES, SCATTER_FORWARD);
};
MOFEM_LOG("FS", Sev::inform) << "Solve projection";
CHKERR solve(D);
auto glob_x = createDMVector(simple->getDM());
auto sub_x = createDMVector(subdm);
auto dummy_dm = create_dummy_dm();
/**
* @brief get TSTheta data operators
*/
auto apply_restrict = [&]() {
auto get_is = [](auto v) {
IS iy;
auto create = [&]() {
int n, ystart;
CHKERR VecGetLocalSize(v, &n);
CHKERR VecGetOwnershipRange(v, &ystart, NULL);
CHKERR ISCreateStride(PETSC_COMM_SELF, n, ystart, 1, &iy);
};
CHK_THROW_MESSAGE(create(), "create is");
return SmartPetscObj<IS>(iy);
};
auto iy = get_is(glob_x);
auto s = createVecScatter(glob_x, PETSC_NULL, glob_x, iy);
DMSubDomainRestrict(simple->getDM(), s, PETSC_NULL, dummy_dm),
"restrict");
Vec X0, Xdot;
CHK_THROW_MESSAGE(DMGetNamedGlobalVector(dummy_dm, "TSTheta_X0", &X0),
"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);
};
CHK_THROW_MESSAGE(forward_ghost(X0), "");
CHK_THROW_MESSAGE(forward_ghost(Xdot), "");
if constexpr (debug) {
MOFEM_LOG("FS", Sev::inform)
<< "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";
CHKERR DMoFEMMeshToLocalVector(simple->getDM(), v, INSERT_VALUES,
SCATTER_REVERSE);
CHKERR solve(sub_x);
for (auto f : {"U", "P", "H", "G", "L"}) {
MOFEM_LOG("WORLD", Sev::verbose) << "Zero field " << f;
CHKERR field_blas->fieldLambdaOnEntities(zero_dofs, f);
}
CHKERR DMoFEMMeshToLocalVector(subdm, sub_x, INSERT_VALUES,
SCATTER_REVERSE);
CHKERR DMoFEMMeshToLocalVector(simple->getDM(), v, INSERT_VALUES,
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);
}
CHKERR DMoFEMMeshToLocalVector(subdm, D, INSERT_VALUES, SCATTER_REVERSE);
}
};
// postprocessing (only for debugging proposes)
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;
post_proc_fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
fe_parent->getOpPtrVector(), {H1});
return fe_parent;
},
QUIET, Sev::noisy);
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",
pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "P",
pip.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "H",
pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "G",
pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
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}},
{},
{}
)
);
auto dm = simple->getDM();
CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
CHKERR post_proc_fe->writeFile("out_projection.h5m");
};
CHKERR solve_projection([](FEMethod *fe_ptr) {
return get_fe_bit(fe_ptr).test(get_current_bit());
});
if constexpr (debug) {
CHKERR post_proc([](FEMethod *fe_ptr) {
return get_fe_bit(fe_ptr).test(get_current_bit());
});
}
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());
}
//! [Data projection]
//! [Push operators to pip]
using UDO = ForcesAndSourcesCore::UserDataOperator;
auto add_parent_field_domain = [&](auto fe, auto op, auto field) {
return setParentDofs(
fe, field, op, bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
QUIET, Sev::noisy);
};
auto add_parent_field_bdy = [&](auto fe, auto op, auto field) {
return setParentDofs(
fe, field, op, bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
QUIET, Sev::noisy);
};
auto test_bit_child = [](FEMethod *fe_ptr) {
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>();
// Push element from reference configuration to current configuration in 3d
// space
auto set_domain_general = [&](auto fe) {
auto &pip = fe->getOpPtrVector();
fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
fe_parent->getOpPtrVector(), {H1});
return fe_parent;
},
QUIET, Sev::noisy);
CHKERR add_parent_field_domain(fe, UDO::OPROW, "U");
pip.push_back(
pip.push_back(new OpCalculateVectorFieldValues<U_FIELD_DIM>("U", u_ptr));
"U", grad_u_ptr));
switch (coord_type) {
case CARTESIAN:
pip.push_back(
"U", div_u_ptr));
break;
pip.push_back(
"U", div_u_ptr));
break;
default:
SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
"Coordinate system not implemented");
}
CHKERR add_parent_field_domain(fe, UDO::OPROW, "H");
pip.push_back(new OpCalculateScalarFieldValuesDot("H", dot_h_ptr));
pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
pip.push_back(
CHKERR add_parent_field_domain(fe, UDO::OPROW, "G");
pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
pip.push_back(
CHKERR add_parent_field_domain(fe, UDO::OPROW, "P");
pip.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
};
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");
pip.push_back(new OpRhsG<false>("G", h_ptr, grad_h_ptr, g_ptr));
CHKERR add_parent_field_domain(fe, UDO::OPROW, "P");
pip.push_back(new OpDomainAssembleScalar(
"P", div_u_ptr, [](const double r, const double, const double) {
return cylindrical(r);
}));
pip.push_back(new OpDomainAssembleScalar(
"P", p_ptr, [](const double r, const double, const double) {
return eps * cylindrical(r);
}));
};
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");
pip.push_back(new OpLhsH_dH<false>("H", u_ptr, h_ptr, grad_g_ptr));
CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G");
pip.push_back(new OpLhsH_dG<false>("H", "G", h_ptr));
}
CHKERR add_parent_field_domain(fe, UDO::OPROW, "G");
{
CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H");
pip.push_back(new OpLhsG_dH<false>("G", "H", h_ptr));
CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G");
pip.push_back(new OpLhsG_dG("G"));
}
CHKERR add_parent_field_domain(fe, UDO::OPROW, "P");
{
CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U");
switch (coord_type) {
case CARTESIAN:
pip.push_back(new OpMixScalarTimesDiv<CARTESIAN>(
"P", "U",
[](const double r, const double, const double) {
return cylindrical(r);
},
true, false));
break;
"P", "U",
[](const double r, const double, const double) {
return cylindrical(r);
},
true, false));
break;
default:
SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
"Coordinate system not implemented");
}
CHKERR add_parent_field_domain(fe, UDO::OPCOL, "P");
pip.push_back(new OpDomainMassP("P", "P", [](double r, double, double) {
return eps * cylindrical(r);
}));
}
};
auto get_block_name = [](auto name) {
return boost::format("%s(.*)") % "WETTING_ANGLE";
};
auto get_blocks = [&](auto &&name) {
return mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
std::regex(name.str()));
};
auto set_boundary_rhs = [&](auto fe) {
auto &pip = fe->getOpPtrVector();
fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
QUIET, Sev::noisy);
CHKERR add_parent_field_bdy(fe, UDO::OPROW, "U");
pip.push_back(new OpCalculateVectorFieldValues<U_FIELD_DIM>("U", u_ptr));
CHKERR add_parent_field_bdy(fe, UDO::OPROW, "L");
pip.push_back(new OpCalculateScalarFieldValues("L", lambda_ptr));
pip.push_back(new OpNormalConstrainRhs("L", u_ptr));
CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpFluidFlux>::add(
pip, mField, "L", {}, "INFLUX",
[](double r, double, double) { return cylindrical(r); }, Sev::inform);
CHKERR add_parent_field_bdy(fe, UDO::OPROW, "U");
pip.push_back(new OpNormalForceRhs("U", lambda_ptr));
auto wetting_block = get_blocks(get_block_name("WETTING_ANGLE"));
if (wetting_block.size()) {
// push operators to the side element which is called from op_bdy_side
auto op_bdy_side =
new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
op_bdy_side->getSideFEPtr()->exeTestHook = test_bit_child;
op_bdy_side->getOpPtrVector(), {H1});
op_bdy_side->getSideFEPtr(), "", UDO::OPSPACE,
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
fe_parent->getOpPtrVector(), {H1});
return fe_parent;
},
QUIET, Sev::noisy);
CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPROW,
"H");
op_bdy_side->getOpPtrVector().push_back(
// push bdy side op
pip.push_back(op_bdy_side);
// push operators for rhs wetting angle
for (auto &b : wetting_block) {
Range force_edges;
std::vector<double> attr_vec;
CHKERR b->getMeshsetIdEntitiesByDimension(
mField.get_moab(), SPACE_DIM - 1, force_edges, true);
b->getAttributes(attr_vec);
if (attr_vec.size() != 1)
SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
"Should be one attribute");
MOFEM_LOG("FS", Sev::inform) << "Wetting angle: " << attr_vec.front();
// need to find the attributes and pass to operator
CHKERR add_parent_field_bdy(fe, UDO::OPROW, "G");
pip.push_back(new OpWettingAngleRhs(
"G", grad_h_ptr, boost::make_shared<Range>(force_edges),
attr_vec.front()));
}
}
};
auto set_boundary_lhs = [&](auto fe) {
auto &pip = fe->getOpPtrVector();
fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
QUIET, Sev::noisy);
CHKERR add_parent_field_bdy(fe, UDO::OPROW, "L");
CHKERR add_parent_field_bdy(fe, UDO::OPCOL, "U");
pip.push_back(new OpNormalConstrainLhs("L", "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>>();
// push operators to the side element which is called from op_bdy_side
auto op_bdy_side =
new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
op_bdy_side->getSideFEPtr()->exeTestHook = test_bit_child;
op_bdy_side->getOpPtrVector(), {H1});
op_bdy_side->getSideFEPtr(), "", UDO::OPSPACE,
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
fe_parent->getOpPtrVector(), {H1});
return fe_parent;
},
QUIET, Sev::noisy);
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(
new OpLoopSideGetDataForSideEle("H", col_ind_ptr, col_diff_base_ptr));
// push bdy side op
pip.push_back(op_bdy_side);
// push operators for lhs wetting angle
for (auto &b : wetting_block) {
Range force_edges;
std::vector<double> attr_vec;
CHKERR b->getMeshsetIdEntitiesByDimension(
mField.get_moab(), SPACE_DIM - 1, force_edges, true);
b->getAttributes(attr_vec);
if (attr_vec.size() != 1)
SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
"Should be one attribute");
MOFEM_LOG("FS", Sev::inform)
<< "wetting angle edges size " << force_edges.size();
CHKERR add_parent_field_bdy(fe, UDO::OPROW, "G");
pip.push_back(new OpWettingAngleLhs(
"G", grad_h_ptr, col_ind_ptr, col_diff_base_ptr,
boost::make_shared<Range>(force_edges), attr_vec.front()));
}
}
};
auto *pip_mng = mField.getInterface<PipelineManager>();
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());
CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule);
CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule);
CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule);
CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule);
pip_mng->getDomainLhsFE()->exeTestHook = test_bit_child;
pip_mng->getDomainRhsFE()->exeTestHook = test_bit_child;
pip_mng->getBoundaryLhsFE()->exeTestHook = test_bit_child;
pip_mng->getBoundaryRhsFE()->exeTestHook = test_bit_child;
}
//! [Push operators to pip]
/**
* @brief Monitor solution
*
* This functions is called by TS kso at the end of each step. It is used
*/
struct Monitor : public FEMethod {
SmartPetscObj<DM> dm, boost::shared_ptr<moab::Core> post_proc_mesh,
boost::shared_ptr<PostProcEleDomainCont> post_proc,
boost::shared_ptr<PostProcEleBdyCont> post_proc_edge,
std::pair<boost::shared_ptr<BoundaryEle>, boost::shared_ptr<VectorDouble>>
p)
: dM(dm), postProcMesh(post_proc_mesh), postProc(post_proc),
postProcEdge(post_proc_edge), liftFE(p.first), liftVec(p.second) {}
MOFEM_LOG("FS", Sev::verbose) << "Monitor";
constexpr int save_every_nth_step = 1;
MOFEM_LOG("FS", Sev::verbose) << "Mesh pre proc";
MoFEM::Interface *m_field_ptr;
auto post_proc_begin =
boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
*m_field_ptr, postProcMesh);
CHKERR DMoFEMPreProcessFiniteElements(dM, post_proc_begin->getFEMethod());
this->getCacheWeakPtr());
this->getCacheWeakPtr());
CHKERR DMoFEMPostProcessFiniteElements(dM, post_proc_end->getFEMethod());
CHKERR post_proc_end->writeFile(
"out_step_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
MOFEM_LOG("FS", Sev::verbose) << "Mesh pre proc done";
}
liftVec->resize(SPACE_DIM, false);
liftVec->clear();
MPI_Allreduce(MPI_IN_PLACE, &(*liftVec)[0], SPACE_DIM, MPI_DOUBLE, MPI_SUM,
MPI_COMM_WORLD);
MOFEM_LOG("FS", Sev::inform)
<< "Step " << ts_step << " time " << ts_t
<< " lift vec x: " << (*liftVec)[0] << " y: " << (*liftVec)[1];
}
private:
boost::shared_ptr<moab::Core> postProcMesh;
boost::shared_ptr<PostProcEleDomainCont> postProc;
boost::shared_ptr<PostProcEleBdyCont> postProcEdge;
boost::shared_ptr<BoundaryEle> liftFE;
boost::shared_ptr<VectorDouble> liftVec;
};
//! [Solve]
using UDO = ForcesAndSourcesCore::UserDataOperator;
auto pip_mng = mField.getInterface<PipelineManager>();
auto create_solver_dm = [&](auto dm) -> SmartPetscObj<DM> {
DM subdm;
auto setup_subdm = [&](auto dm) {
auto bit_mng = mField.getInterface<BitRefManager>();
auto dm = simple->getDM();
auto level_ents_ptr = boost::make_shared<Range>();
CHKERR bit_mng->getEntitiesByRefLevel(
bit(get_current_bit()), BitRefLevel().set(), *level_ents_ptr);
CHKERR DMCreate(mField.get_comm(), &subdm);
CHKERR DMSetType(subdm, "DMMOFEM");
CHKERR DMMoFEMCreateSubDM(subdm, dm, "SUB_SOLVER");
CHKERR DMMoFEMAddElement(subdm, simple->getDomainFEName());
CHKERR DMMoFEMAddElement(subdm, simple->getBoundaryFEName());
CHKERR DMMoFEMSetSquareProblem(subdm, PETSC_TRUE);
CHKERR DMMoFEMSetDestroyProblem(subdm, PETSC_FALSE);
for (auto f : {"U", "P", "H", "G", "L"}) {
CHKERR DMMoFEMAddSubFieldRow(subdm, f, level_ents_ptr);
CHKERR DMMoFEMAddSubFieldCol(subdm, f, level_ents_ptr);
}
CHKERR DMSetUp(subdm);
};
CHK_THROW_MESSAGE(setup_subdm(dm), "create subdm");
return SmartPetscObj<DM>(subdm);
};
auto get_fe_post_proc = [&](auto post_proc_mesh) {
auto add_parent_field_domain = [&](auto fe, auto op, auto field) {
return setParentDofs(
fe, field, op, bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
return fe_parent;
},
QUIET, Sev::noisy);
};
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});
post_proc_fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
fe_parent->getOpPtrVector(), {H1});
return fe_parent;
},
QUIET, Sev::noisy);
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(
new OpCalculateScalarFieldValues("H", h_ptr));
post_proc_fe->getOpPtrVector().push_back(
CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "P");
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("P", p_ptr));
CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "G");
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("G", g_ptr));
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(
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(
fe, field, op, bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
new BoundaryParentEle(mField));
return fe_parent;
},
QUIET, Sev::noisy);
};
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(
};
CHKERR setParentDofs(
post_proc_fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
new BoundaryParentEle(mField));
return fe_parent;
},
QUIET, Sev::noisy);
struct OpGetNormal : public BoundaryEleOp {
OpGetNormal(boost::shared_ptr<VectorDouble> l_ptr,
boost::shared_ptr<MatrixDouble> n_ptr)
: BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE), ptrL(l_ptr),
ptrNormal(n_ptr) {}
auto t_l = getFTensor0FromVec(*ptrL);
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(
new OpCalculateScalarFieldValues("L", lambda_ptr));
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(
new OpCalculateScalarFieldValues("P", p_ptr));
op_ptr->doWorkRhsHook = [&](DataOperator *base_op_ptr, int side,
auto op_ptr = static_cast<BoundaryEleOp *>(base_op_ptr);
};
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(
fe, field, op, bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
new BoundaryParentEle(mField));
return fe_parent;
},
QUIET, Sev::noisy);
};
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>();
CHKERR setParentDofs(
fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
new BoundaryParentEle(mField));
return fe_parent;
},
QUIET, Sev::noisy);
std::vector<const CubitMeshSets *> vec_ptr;
CHKERR mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
std::regex("LIFT"), vec_ptr);
for (auto m_ptr : vec_ptr) {
auto meshset = m_ptr->getMeshset();
Range ents;
CHKERR mField.get_moab().get_entities_by_dimension(meshset, SPACE_DIM - 1,
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(
new OpCalculateScalarFieldValues("L", p_ptr));
fe->getOpPtrVector().push_back(
new OpCalculateLift("L", p_ptr, lift_ptr, ents_ptr));
return std::make_pair(fe, lift_ptr);
};
auto set_post_proc_monitor = [&](auto ts) {
DM dm;
CHKERR TSGetDM(ts, &dm);
boost::shared_ptr<FEMethod> null_fe;
auto post_proc_mesh = boost::make_shared<moab::Core>();
auto monitor_ptr = boost::make_shared<Monitor>(
SmartPetscObj<DM>(dm, true), post_proc_mesh,
get_fe_post_proc(post_proc_mesh), get_bdy_post_proc_fe(post_proc_mesh),
get_lift_fe());
CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
null_fe, monitor_ptr);
};
auto dm = simple->getDM();
auto ts = createTS(mField.get_comm());
CHKERR TSSetDM(ts, dm);
auto ts_pre_post_proc = boost::make_shared<TSPrePostProc>();
tsPrePostProc = ts_pre_post_proc;
if (auto ptr = tsPrePostProc.lock()) {
ptr->fsRawPtr = this;
ptr->solverSubDM = create_solver_dm(simple->getDM());
ptr->globSol = createDMVector(dm);
CHKERR DMoFEMMeshToLocalVector(dm, ptr->globSol, INSERT_VALUES,
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);
// Add monitor to time solver
double ftime = 1;
CHKERR TSSetFromOptions(ts);
CHKERR ptr->tsSetUp(ts);
CHKERR TSSetUp(ts);
auto print_fields_in_section = [&]() {
auto section = mField.getInterface<ISManager>()->sectionCreate(
simple->getProblemName());
PetscInt num_fields;
CHKERR PetscSectionGetNumFields(section, &num_fields);
for (int f = 0; f < num_fields; ++f) {
const char *field_name;
CHKERR PetscSectionGetFieldName(section, f, &field_name);
MOFEM_LOG("FS", Sev::inform)
<< "Field " << f << " " << std::string(field_name);
}
};
CHKERR print_fields_in_section();
CHKERR TSSolve(ts, ptr->globSol);
}
}
int main(int argc, char *argv[]) {
#ifdef PYTHON_INIT_SURFACE
Py_Initialize();
#endif
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
// Add logging channel for example
auto core_log = logging::core::get();
core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "FS"));
LogManager::setLog("FS");
MOFEM_LOG_TAG("FS", "free surface");
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
//! [Register MoFEM discrete manager in PETSc
//! [Create MoAB]
moab::Core mb_instance; ///< mesh database
moab::Interface &moab = mb_instance; ///< mesh database interface
//! [Create MoAB]
//! [Create MoFEM]
MoFEM::Core core(moab); ///< finite element database
MoFEM::Interface &m_field = core; ///< finite element database interface
//! [Create MoFEM]
//! [FreeSurface]
FreeSurface ex(m_field);
CHKERR ex.runProblem();
//! [FreeSurface]
}
#ifdef PYTHON_INIT_SURFACE
if (Py_FinalizeEx() < 0) {
exit(120);
}
#endif
}
std::vector<Range>
auto &moab = mField.get_moab();
auto bit_mng = mField.getInterface<BitRefManager>();
auto comm_mng = mField.getInterface<CommInterface>();
Range vertices;
CHK_THROW_MESSAGE(bit_mng->getEntitiesByTypeAndRefLevel(
bit(0), BitRefLevel().set(), MBVERTEX, vertices),
"can not get vertices on bit 0");
auto &dofs_mi = mField.get_dofs()->get<Unique_mi_tag>();
auto field_bit_number = mField.get_field_bit_number("H");
Range plus_range, minus_range;
std::vector<EntityHandle> plus, minus;
// get vertices on level 0 on plus and minus side
for (auto p = vertices.pair_begin(); p != vertices.pair_end(); ++p) {
const auto f = p->first;
const auto s = p->second;
// Lowest Dof UId for given field (field bit number) on entity f
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();
if (v > 0)
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());
}
MOFEM_TAG_AND_LOG("SYNC", Sev::noisy, "FS")
<< "Plus range " << plus_range << endl;
MOFEM_TAG_AND_LOG("SYNC", Sev::noisy, "FS")
<< "Minus range " << minus_range << endl;
auto get_elems = [&](auto &ents, auto bit, auto mask) {
Range adj;
CHK_MOAB_THROW(moab.get_adjacencies(ents, SPACE_DIM, false, adj,
moab::Interface::UNION),
"can not get adjacencies");
CHK_THROW_MESSAGE(bit_mng->filterEntitiesByRefLevel(bit, mask, adj),
"can not filter elements with bit 0");
return adj;
};
CHKERR comm_mng->synchroniseEntities(plus_range);
CHKERR comm_mng->synchroniseEntities(minus_range);
std::vector<Range> ele_plus(nb_levels), ele_minus(nb_levels);
ele_plus[0] = get_elems(plus_range, bit(0), BitRefLevel().set());
ele_minus[0] = get_elems(minus_range, bit(0), BitRefLevel().set());
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);
c = c.subset_by_dimension(SPACE_DIM);
};
for (auto l = 1; l != nb_levels; ++l) {
CHK_THROW_MESSAGE(get_children(ele_plus[l - 1], ele_plus[l]),
"get children");
CHK_THROW_MESSAGE(get_children(ele_minus[l - 1], ele_minus[l]),
"get children");
}
auto get_level = [&](auto &p, auto &m, auto z, auto bit, auto mask) {
bit_mng->getEntitiesByDimAndRefLevel(bit, mask, SPACE_DIM, l),
"can not get vertices on bit");
l = subtract(l, p);
l = subtract(l, m);
for (auto f = 0; f != z; ++f) {
Range conn;
CHK_MOAB_THROW(moab.get_connectivity(l, conn, true), "");
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(conn);
l = get_elems(conn, bit, mask);
}
return l;
};
std::vector<Range> vec_levels(nb_levels);
for (auto l = nb_levels - 1; l >= 0; --l) {
vec_levels[l] = get_level(ele_plus[l], ele_minus[l], 2 * overlap, bit(l),
BitRefLevel().set());
}
if constexpr (debug) {
if (mField.get_comm_size() == 1) {
for (auto l = 0; l != nb_levels; ++l) {
std::string name = (boost::format("out_r%d.h5m") % l).str();
CHK_THROW_MESSAGE(save_range(mField.get_moab(), name, vec_levels[l]),
"save mesh");
}
}
}
return vec_levels;
}
auto bit_mng = mField.getInterface<BitRefManager>();
auto prb_mng = mField.getInterface<ProblemsManager>();
BitRefLevel start_mask;
for (auto s = 0; s != get_start_bit(); ++s)
start_mask[s] = true;
// store prev_level
Range prev_level;
CHKERR bit_mng->getEntitiesByRefLevel(bit(get_current_bit()),
BitRefLevel().set(), prev_level);
Range prev_level_skin;
CHKERR bit_mng->getEntitiesByRefLevel(bit(get_skin_parent_bit()),
BitRefLevel().set(), prev_level_skin);
// reset bit ref levels
CHKERR bit_mng->lambdaBitRefLevel(
[&](EntityHandle ent, BitRefLevel &bit) { bit &= start_mask; });
CHKERR bit_mng->setNthBitRefLevel(prev_level, get_projection_bit(), true);
CHKERR bit_mng->setNthBitRefLevel(prev_level_skin, get_skin_projection_bit(),
true);
// set refinement levels
auto set_levels = [&](auto &&
vec_levels /*entities are refined on each level*/) {
// start with zero level, which is the coarsest mesh
Range level0;
CHKERR bit_mng->getEntitiesByRefLevel(bit(0), BitRefLevel().set(), level0);
CHKERR bit_mng->setNthBitRefLevel(level0, get_start_bit(), true);
// get lower dimension entities
auto get_adj = [&](auto ents) {
Range conn;
CHK_MOAB_THROW(mField.get_moab().get_connectivity(ents, conn, true),
"get conn");
for (auto d = 1; d != SPACE_DIM; ++d) {
CHK_MOAB_THROW(mField.get_moab().get_adjacencies(
ents.subset_by_dimension(SPACE_DIM), d, false, conn,
moab::Interface::UNION),
"get adj");
}
ents.merge(conn);
return ents;
};
// set bit levels
for (auto l = 1; l != nb_levels; ++l) {
Range level_prev;
CHKERR bit_mng->getEntitiesByDimAndRefLevel(bit(get_start_bit() + l - 1),
BitRefLevel().set(),
SPACE_DIM, level_prev);
Range parents;
CHKERR bit_mng->updateRangeByParent(vec_levels[l], parents);
// subtract entities from previous level, which are refined, so should be
// not there
level_prev = subtract(level_prev, parents);
// and instead add their children
level_prev.merge(vec_levels[l]);
// set bit to each level
CHKERR bit_mng->setNthBitRefLevel(level_prev, get_start_bit() + l, true);
}
// set bit levels to lower dimension entities
for (auto l = 1; l != nb_levels; ++l) {
Range level;
CHKERR bit_mng->getEntitiesByDimAndRefLevel(
bit(get_start_bit() + l), BitRefLevel().set(), SPACE_DIM, level);
level = get_adj(level);
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(level);
CHKERR bit_mng->setNthBitRefLevel(level, get_start_bit() + l, true);
}
};
// resolve skin between refined levels
auto set_skins = [&]() {
moab::Skinner skinner(&mField.get_moab());
ParallelComm *pcomm =
ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
// get skin of bit level
auto get_bit_skin = [&](BitRefLevel bit, BitRefLevel mask) {
Range bit_ents;
bit, mask, SPACE_DIM, bit_ents),
"can't get bit level");
Range bit_skin;
CHK_MOAB_THROW(skinner.find_skin(0, bit_ents, false, bit_skin),
"can't get skin");
return bit_skin;
};
auto get_level_skin = [&]() {
Range skin;
BitRefLevel bit_prev;
for (auto l = 1; l != nb_levels; ++l) {
auto skin_level_mesh = get_bit_skin(bit(l), BitRefLevel().set());
// filter (remove) all entities which are on partitions boarders
CHKERR pcomm->filter_pstatus(skin_level_mesh,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, nullptr);
auto skin_level =
get_bit_skin(bit(get_start_bit() + l), BitRefLevel().set());
skin_level = subtract(skin_level,
skin_level_mesh); // get only internal skins, not
// on the body boundary
// get lower dimension adjacencies (FIXME: add edges if 3D)
Range skin_level_verts;
CHKERR mField.get_moab().get_connectivity(skin_level, skin_level_verts,
true);
skin_level.merge(skin_level_verts);
// remove previous level
bit_prev.set(l - 1);
Range level_prev;
CHKERR bit_mng->getEntitiesByRefLevel(bit_prev, BitRefLevel().set(),
level_prev);
skin.merge(subtract(skin_level, level_prev));
}
return skin;
};
auto resolve_shared = [&](auto &&skin) {
Range tmp_skin = skin;
map<int, Range> map_procs;
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
tmp_skin, &map_procs);
Range from_other_procs; // entities which also exist on other processors
for (auto &m : map_procs) {
if (m.first != mField.get_comm_rank()) {
from_other_procs.merge(m.second);
}
}
auto common = intersect(
skin, from_other_procs); // entities which are on internal skin
skin.merge(from_other_procs);
// entities which are on processors boards, and several processors are not
// true skin.
if (!common.empty()) {
// skin is internal exist on other procs
skin = subtract(skin, common);
}
return skin;
};
// get parents of entities
auto get_parent_level_skin = [&](auto skin) {
Range skin_parents;
CHKERR bit_mng->updateRangeByParent(
skin.subset_by_dimension(SPACE_DIM - 1), skin_parents);
Range skin_parent_verts;
CHKERR mField.get_moab().get_connectivity(skin_parents, skin_parent_verts,
true);
skin_parents.merge(skin_parent_verts);
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
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);
CHKERR bit_mng->setNthBitRefLevel(child_skin, get_skin_child_bit(), true);
CHKERR bit_mng->setNthBitRefLevel(parent_skin, get_skin_parent_bit(), true);
};
// take last level, remove childs on boarder, and set bit
auto set_current = [&]() {
Range last_level;
CHKERR bit_mng->getEntitiesByRefLevel(bit(get_start_bit() + nb_levels - 1),
BitRefLevel().set(), last_level);
Range skin_child;
CHKERR bit_mng->getEntitiesByRefLevel(bit(get_skin_child_bit()),
BitRefLevel().set(), skin_child);
last_level = subtract(last_level, skin_child);
CHKERR bit_mng->setNthBitRefLevel(last_level, get_current_bit(), true);
};
// set bits to levels
// set bits to skin
CHKERR set_skins();
// set current level bit
CHKERR set_current();
if constexpr (debug) {
if (mField.get_comm_size() == 1) {
for (auto l = 0; l != nb_levels; ++l) {
std::string name = (boost::format("out_level%d.h5m") % l).str();
CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_start_bit() + l),
BitRefLevel().set(), name.c_str(), "MOAB",
"PARALLEL=WRITE_PART");
}
CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_current_bit()),
BitRefLevel().set(), "current_bit.h5m",
"MOAB", "PARALLEL=WRITE_PART");
CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_projection_bit()),
BitRefLevel().set(), "projection_bit.h5m",
"MOAB", "PARALLEL=WRITE_PART");
CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_skin_child_bit()),
BitRefLevel().set(), "skin_child_bit.h5m",
"MOAB", "PARALLEL=WRITE_PART");
CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_skin_parent_bit()),
BitRefLevel().set(), "skin_parent_bit.h5m",
"MOAB", "PARALLEL=WRITE_PART");
}
}
};
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) {
/**
* @brief Collect data from parent elements to child
*/
boost::function<void(boost::shared_ptr<ForcesAndSourcesCore>, int)>
add_parent_level =
[&](boost::shared_ptr<ForcesAndSourcesCore> parent_fe_pt, int level) {
// Evaluate if not last parent element
if (level > 0) {
// Create domain parent FE
auto fe_ptr_current = get_elem();
// Call next level
add_parent_level(
boost::dynamic_pointer_cast<ForcesAndSourcesCore>(
fe_ptr_current),
level - 1);
// Add data to curent fe level
if (op == ForcesAndSourcesCore::UserDataOperator::OPSPACE) {
// Only base
parent_fe_pt->getOpPtrVector().push_back(
H1, op, fe_ptr_current,
BitRefLevel().set(), BitRefLevel().set(),
child_ent_bit, BitRefLevel().set(),
verbosity, sev));
} else {
// Filed data
parent_fe_pt->getOpPtrVector().push_back(
field_name, op, fe_ptr_current,
BitRefLevel().set(), BitRefLevel().set(),
child_ent_bit, BitRefLevel().set(),
verbosity, sev));
}
}
};
add_parent_level(boost::dynamic_pointer_cast<ForcesAndSourcesCore>(fe_top),
}
if (auto ptr = tsPrePostProc.lock()) {
/**
* @brief cut-off values at nodes, i.e. abs("H") <= 1
*
*/
auto cut_off_dofs = [&]() {
auto &m_field = ptr->fsRawPtr->mField;
Range current_verts;
auto bit_mng = m_field.getInterface<BitRefManager>();
CHKERR bit_mng->getEntitiesByTypeAndRefLevel(
bit(get_current_bit()), BitRefLevel().set(), MBVERTEX, current_verts);
auto cut_off_verts = [&](boost::shared_ptr<FieldEntity> ent_ptr) {
for (auto &h : ent_ptr->getEntFieldData()) {
h = cut_off(h);
}
};
auto field_blas = m_field.getInterface<FieldBlas>();
CHKERR field_blas->fieldLambdaOnEntities(cut_off_verts, "H",
&current_verts);
};
CHKERR cut_off_dofs();
}
if (auto ptr = tsPrePostProc.lock()) {
MOFEM_LOG("FS", Sev::inform) << "Run step pre proc";
auto &m_field = ptr->fsRawPtr->mField;
auto simple = m_field.getInterface<Simple>();
auto bit_mng = m_field.getInterface<BitRefManager>();
auto bc_mng = m_field.getInterface<BcManager>();
auto field_blas = m_field.getInterface<FieldBlas>();
auto opt = m_field.getInterface<OperatorsTester>();
auto pip_mng = m_field.getInterface<PipelineManager>();
auto prb_mng = m_field.getInterface<ProblemsManager>();
// get vector norm
auto get_norm = [&](auto x) {
double nrm;
CHKERR VecNorm(x, NORM_2, &nrm);
return nrm;
};
// refine problem and project data, including theta data
auto refine_problem = [&]() {
MOFEM_LOG("FS", Sev::inform) << "Refine problem";
CHKERR ptr->fsRawPtr->refineMesh(refine_overlap);
CHKERR ptr->fsRawPtr->projectData();
};
// set new jacobin operator, since problem and thus tangent matrix size has
// changed
auto set_jacobian_operators = [&]() {
ptr->subB = createDMMatrix(ptr->solverSubDM);
CHKERR KSPReset(ptr->subKSP);
};
// set new solution
auto set_solution = [&]() {
MOFEM_LOG("FS", Sev::inform) << "Set solution";
PetscObjectState state;
// Record the state, and set it again. This is to fool PETSc that solution
// vector is not updated. Otherwise PETSc will treat every step as a first
// step.
// globSol is updated as result mesh refinement - this is not really set
// a new solution.
CHKERR PetscObjectStateGet(getPetscObject(ptr->globSol.get()), &state);
CHKERR DMoFEMMeshToLocalVector(simple->getDM(), ptr->globSol,
INSERT_VALUES, SCATTER_FORWARD);
CHKERR PetscObjectStateSet(getPetscObject(ptr->globSol.get()), state);
MOFEM_LOG("FS", Sev::verbose)
<< "Set solution, vector norm " << get_norm(ptr->globSol);
};
PetscBool is_theta;
PetscObjectTypeCompare((PetscObject)ts, TSTHETA, &is_theta);
if (is_theta) {
CHKERR refine_problem(); // refine problem
CHKERR set_jacobian_operators(); // set new jacobian
CHKERR set_solution(); // set solution
} else {
SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
"Sorry, only TSTheta handling is implemented");
}
// Need barriers, somme functions in TS solver need are called collectively
// and requite the same state of variables
PetscBarrier((PetscObject)ts);
MOFEM_TAG_AND_LOG("SYNC", Sev::verbose, "FS") << "PreProc done";
MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::verbose);
}
}
if (auto ptr = tsPrePostProc.lock()) {
auto &m_field = ptr->fsRawPtr->mField;
MOFEM_TAG_AND_LOG("SYNC", Sev::verbose, "FS") << "PostProc done";
MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::verbose);
}
return 0;
}
MoFEMErrorCode TSPrePostProc::tsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t,
Vec f, void *ctx) {
if (auto ptr = tsPrePostProc.lock()) {
auto sub_u = ptr->getSubVector();
auto sub_u_t = vectorDuplicate(sub_u);
auto sub_f = vectorDuplicate(sub_u);
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 TsSetIFunction(ts, t, sub_u, sub_u_t, sub_f, ptr->tsCtxPtr.get());
CHKERR VecScatterBegin(scatter, sub_f, f, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecScatterEnd(scatter, sub_f, f, INSERT_VALUES, SCATTER_FORWARD);
}
}
MoFEMErrorCode TSPrePostProc::tsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t,
PetscReal a, Mat A, Mat B,
void *ctx) {
if (auto ptr = tsPrePostProc.lock()) {
auto sub_u = ptr->getSubVector();
auto sub_u_t = vectorDuplicate(sub_u);
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 TsSetIJacobian(ts, t, sub_u, sub_u_t, a, ptr->subB, ptr->subB,
ptr->tsCtxPtr.get());
}
}
MoFEMErrorCode TSPrePostProc::tsMonitor(TS ts, PetscInt step, PetscReal t,
Vec u, void *ctx) {
if (auto ptr = tsPrePostProc.lock()) {
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);
MOFEM_LOG("FS", Sev::verbose)
<< "u norm " << get_norm(u) << " u sub nom " << get_norm(sub_u);
CHKERR TsMonitorSet(ts, step, t, sub_u, ptr->tsCtxPtr.get());
}
}
if (auto ptr = tsPrePostProc.lock()) {
MOFEM_LOG("FS", Sev::verbose) << "SetUP sub PC";
ptr->subKSP = createKSP(ptr->fsRawPtr->mField.get_comm());
CHKERR KSPSetFromOptions(ptr->subKSP);
}
};
MoFEMErrorCode TSPrePostProc::pcApply(PC pc, Vec pc_f, Vec pc_x) {
if (auto ptr = tsPrePostProc.lock()) {
auto sub_x = ptr->getSubVector();
auto sub_f = vectorDuplicate(sub_x);
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);
}
};
if (auto ptr = tsPrePostProc.lock()) {
auto prb_ptr = ptr->fsRawPtr->mField.get_problem("SUB_SOLVER");
if (auto sub_data = prb_ptr->getSubData()) {
auto is = sub_data->getSmartColIs();
VecScatter s;
if (fr == R) {
CHK_THROW_MESSAGE(VecScatterCreate(x, PETSC_NULL, y, is, &s),
"crate scatter");
} else {
CHK_THROW_MESSAGE(VecScatterCreate(x, is, y, PETSC_NULL, &s),
"crate scatter");
}
}
}
}
}
auto &m_field = fsRawPtr->mField;
auto simple = m_field.getInterface<Simple>();
auto dm = simple->getDM();
CHKERR TSSetIFunction(ts, globRes, tsSetIFunction, nullptr);
CHKERR TSSetIJacobian(ts, PETSC_NULL, PETSC_NULL, tsSetIJacobian, nullptr);
CHKERR TSMonitorSet(ts, tsMonitor, fsRawPtr, PETSC_NULL);
SNES snes;
CHKERR TSGetSNES(ts, &snes);
auto snes_ctx_ptr = getDMSnesCtx(dm);
auto set_section_monitor = [&](auto snes) {
CHKERR SNESMonitorSet(snes,
(MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
(void *)(snes_ctx_ptr.get()), nullptr);
};
CHKERR set_section_monitor(snes);
auto ksp = createKSP(m_field.get_comm());
CHKERR KSPSetType(ksp, KSPPREONLY); // Run KSP internally in ShellPC
auto sub_pc = createPC(fsRawPtr->mField.get_comm());
CHKERR PCSetType(sub_pc, PCSHELL);
CHKERR PCShellSetSetUp(sub_pc, pcSetup);
CHKERR PCShellSetApply(sub_pc, pcApply);
CHKERR KSPSetPC(ksp, sub_pc);
CHKERR SNESSetKSP(snes, ksp);
CHKERR TSSetPreStep(ts, tsPreProc);
CHKERR TSSetPostStep(ts, tsPostProc);
}
static Index< 'p', 3 > p
std::string param_file
#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)
Definition acoustic.cpp:69
static char help[]
int main()
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
PostProcEleByDim< SPACE_DIM >::SideEle SideEle
constexpr double a
static const double eps
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
BoundaryEle::UserDataOperator BoundaryEleOp
Kronecker Delta class symmetric.
@ QUIET
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
#define BITREFLEVEL_SIZE
max number of refinements
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
@ H1
continuous field
Definition definitions.h:85
@ NOSPACE
Definition definitions.h:83
#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
Definition definitions.h:34
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_INVALID_DATA
Definition definitions.h:36
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
@ LAST_COORDINATE_SYSTEM
@ CYLINDRICAL
@ CARTESIAN
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
constexpr int BASE_DIM
constexpr int order
static const bool debug
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
auto cylindrical
double kappa
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, U_FIELD_DIM > OpDomainMassU
auto bubble_device
double mu_diff
auto save_range
auto get_M2
constexpr int U_FIELD_DIM
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
static boost::weak_ptr< TSPrePostProc > tsPrePostProc
auto get_M2_dh
FormsIntegrators< BoundaryEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, 1 > OpBoundaryMassL
@ R
@ F
auto integration_rule
auto get_M_dh
auto kernel_eye
double mu_m
double lambda
surface tension
PostProcBrokenMeshInMoabBaseCont< DomainEle > PostProcEleDomainCont
OpDomainMassH OpDomainMassG
auto get_M3_dh
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesVector< BASE_DIM, SPACE_DIM, 1 > OpDomainAssembleVector
double rho_diff
auto init_h
Initialisation function.
constexpr auto t_kd
auto d_phase_function_h
auto get_fe_bit
constexpr int SPACE_DIM
double rho_ave
FTensor::Index< 'l', SPACE_DIM > l
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpDomainAssembleScalar
int coord_type
double W
auto kernel_oscillation
PostProcBrokenMeshInMoabBaseCont< BoundaryEle > PostProcEleBdyCont
double eta2
auto get_skin_parent_bit
auto get_M0
constexpr bool debug
auto get_f_dh
auto get_M3
double rho_p
double mu_p
auto get_start_bit
OpDomainMassH OpDomainMassP
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMixScalarTimesDiv< SPACE_DIM, COORD_TYPE > OpMixScalarTimesDiv
int nb_levels
auto wetting_angle
double h
auto get_M
double mu_ave
double eps
auto get_dofs_ents_all
get entities of dofs in the problem - used for debugging
auto phase_function
auto get_skin_projection_bit
auto get_global_size
auto get_dofs_ents_by_field_name
get entities of dofs in the problem - used for debugging
ElementsAndOps< SPACE_DIM >::DomianParentEle DomianParentEle
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, 1 > OpDomainMassH
auto get_D
auto wetting_angle_sub_stepping
int refine_overlap
auto get_M0_dh
auto my_max
auto cut_off
auto get_skin_child_bit
int order
approximation order
SideEle::UserDataOperator SideOp
double a0
double rho_m
auto get_current_bit
dofs bit used to do calculations
double eta
auto bit
double md
auto capillary_tube
FormsIntegrators< BoundaryEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpBoundaryAssembleScalar
auto get_f
BoundaryNaturalBC::OpFlux< NaturalMeshsetType< BLOCKSET >, 1, 1 > OpFluidFlux
auto get_projection_bit
auto d_cut_off
auto my_min
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition DMMoFEM.cpp:219
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:483
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition DMMoFEM.cpp:442
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.
Definition DMMoFEM.cpp:118
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:542
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition DMMoFEM.cpp:509
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition DMMoFEM.cpp:242
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:47
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:572
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1003
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
Definition DMMoFEM.cpp:400
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition DMMoFEM.cpp:275
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition DMMoFEM.hpp:988
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:532
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
IntegrationType
Form integrator integration types.
AssemblyType
[Storage and set boundary conditions]
#define MOFEM_LOG(channel, severity)
Log.
SeverityLevel
Severity levels.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
MoFEMErrorCode getMeshset(const int ms_id, const unsigned int cubit_bc_type, EntityHandle &meshset) const
get meshset from CUBIT Id and CUBIT type
MoFEMErrorCode buildSubProblem(const std::string out_name, const std::vector< std::string > &fields_row, const std::vector< std::string > &fields_col, const std::string main_problem, const bool square_matrix=true, const map< std::string, boost::shared_ptr< Range > > *entityMapRow=nullptr, const map< std::string, boost::shared_ptr< Range > > *entityMapCol=nullptr, int verb=VERBOSE)
build sub problem
auto marker
set bit to marker
auto test_bit_child
lambda function used to select elements on which finite element pipelines are executed.
auto bit
set bit
constexpr double a0
FTensor::Index< 'i', SPACE_DIM > i
constexpr CoordinateTypes coord_type
static double lambda
const double c
speed of light (cm/ns)
double D
const double v
phase velocity of light in medium (cm/ns)
constexpr int nb_levels
Definition level_set.cpp:58
PipelineManager::ElementsAndOpsByDim< FE_DIM >::DomianParentEle DomianParentEle
Definition level_set.cpp:38
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
double tol
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
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.
Definition DMMoFEM.cpp:1042
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobian in TS solver.
Definition TsCtx.cpp:165
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition TsCtx.cpp:259
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition DMMoFEM.hpp:1045
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition DMMoFEM.cpp:424
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
Definition SnesCtx.cpp:232
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.
Definition TsCtx.cpp:56
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.
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition DMMoFEM.hpp:1031
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
auto getProblemPtr(DM dm)
get problem pointer from DM
Definition DMMoFEM.hpp:976
double w(const double g, const double t)
float d
Definition sdf_hertz.py:5
int r
Definition sdf.py:8
constexpr IntegrationType I
constexpr AssemblyType A
double h
int save_every_nth_step
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr double t
plate stiffness
Definition plate.cpp:59
constexpr auto field_name
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition radiation.cpp:29
constexpr double g
MoFEM::Interface & mField
Definition plastic.cpp:235
Range findParentsToRefine(Range ents, BitRefLevel level, BitRefLevel mask)
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode runProblem()
[Run programme]
FreeSurface(MoFEM::Interface &m_field)
MoFEMErrorCode projectData()
[Boundary condition]
MoFEMErrorCode boundaryCondition()
[Set up problem]
MoFEMErrorCode readMesh()
[Run programme]
MoFEMErrorCode setParentDofs(boost::shared_ptr< FEMethod > fe_top, std::string field_name, ForcesAndSourcesCore::UserDataOperator::OpType op, BitRefLevel child_ent_bit, boost::function< boost::shared_ptr< ForcesAndSourcesCore >()> get_elem, int verbosity, LogManager::SeverityLevel sev)
Create hierarchy of elements run on parents levels.
std::vector< Range > findEntitiesCrossedByPhaseInterface(size_t overlap)
Find entities on refinement levels.
MoFEMErrorCode assembleSystem()
[Data projection]
MoFEM::Interface & mField
MoFEMErrorCode refineMesh(size_t overlap)
MoFEMErrorCode makeRefProblem()
MoFEMErrorCode solveSystem()
[Solve]
Add operators pushing bases from local to physical configuration.
boost::weak_ptr< CacheTuple > getCacheWeakPtr() const
Get the cache weak ptr object.
Simple interface for fast problem set-up.
Definition BcManager.hpp:25
Managing BitRefLevels.
Managing BitRefLevels.
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
Core (interface) class.
Definition Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:112
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)
Basic algebra on fields.
Definition FieldBlas.hpp:21
@ 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.
Definition ISManager.hpp:23
Interface for managing meshsets containing materials and boundary conditions.
Assembly methods.
Definition Natural.hpp:65
Operator to project base functions from parent entity to child.
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
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.
Definition Simple.hpp:27
intrusive_ptr for managing petsc objects
PetscReal ts_t
time
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
Monitor(SmartPetscObj< DM > dm, boost::shared_ptr< PostProcFace > post_proc, boost::shared_ptr< FEMethod > reaction_fe, std::vector< boost::shared_ptr< ScalingMethod > > smv)
SmartPetscObj< DM > dM
MoFEMErrorCode postProcess()
function is run at the end of loop
boost::shared_ptr< VectorDouble > liftVec
boost::shared_ptr< BoundaryEle > liftFE
boost::shared_ptr< PostProcFace > 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
Example * fsRawPtr
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.
TSPrePostProc()=default
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
auto save_range
NaturalBC< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS > BoundaryNaturalBC
[Body and heat source]