v0.14.0
VEC-5: Phase field model and Navier-Stokes

Implentaton of Navier-Stokes-Cahn-Hilliard model is based on [48].

Note
Prerequisites of this tutorial include VEC-0: Linear elasticity


Note
Intended learning outcome:
  • handling multiple fields
  • setting inital boundary conditions
  • material time derivative
  • hyperbolic equation

Example solution

Plain program

/**
* \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";
#undef PYTHON_INIT_SURFACE
#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) {
MoFEMFunctionBegin;
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();
CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
}
MoFEMFunctionReturn(0);
};
MoFEMErrorCode evalSurface(
double x, double y, double z, double eta, double &s
) {
MoFEMFunctionBegin;
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();
CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
}
MoFEMFunctionReturn(0);
}
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 ElementsAndOps = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>;
using DomainEle = ElementsAndOps<SPACE_DIM>::DomainEle;
using DomianParentEle = ElementsAndOps<SPACE_DIM>::DomianParentEle;
using DomainEleOp = DomainEle::UserDataOperator;
using BoundaryEle = ElementsAndOps<SPACE_DIM>::BoundaryEle;
using BoundaryParentEle = ElementsAndOps<SPACE_DIM>::BoundaryParentEle;
using BoundaryEleOp = BoundaryEle::UserDataOperator;
using SideEle = ElementsAndOps<SPACE_DIM>::FaceSideEle;
using SideOp = SideEle::UserDataOperator;
using PostProcEleDomain = PostProcBrokenMeshInMoab<DomainEle>;
using PostProcEleDomainCont = PostProcBrokenMeshInMoabBaseCont<DomainEle>;
using PostProcEleBdyCont = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
using EntData = EntitiesFieldData::EntData;
using AssemblyDomainEleOp = FormsIntegrators<DomainEleOp>::Assembly<A>::OpBase;
using AssemblyBoundaryEleOp =
FormsIntegrators<BoundaryEleOp>::Assembly<A>::OpBase;
using OpDomainMassU = FormsIntegrators<DomainEleOp>::Assembly<A>::BiLinearForm<
I>::OpMass<BASE_DIM, U_FIELD_DIM>;
using OpDomainMassH = FormsIntegrators<DomainEleOp>::Assembly<A>::BiLinearForm<
I>::OpMass<BASE_DIM, 1>;
using OpDomainMassP = OpDomainMassH;
using OpDomainMassG = OpDomainMassH;
using OpBoundaryMassL = FormsIntegrators<BoundaryEleOp>::Assembly<
A>::BiLinearForm<I>::OpMass<BASE_DIM, 1>;
using OpDomainAssembleVector = FormsIntegrators<DomainEleOp>::Assembly<
A>::LinearForm<I>::OpBaseTimesVector<BASE_DIM, SPACE_DIM, 1>;
using OpDomainAssembleScalar = FormsIntegrators<DomainEleOp>::Assembly<
A>::LinearForm<I>::OpBaseTimesScalar<BASE_DIM>;
using OpBoundaryAssembleScalar = FormsIntegrators<BoundaryEleOp>::Assembly<
A>::LinearForm<I>::OpBaseTimesScalar<BASE_DIM>;
template <CoordinateTypes COORD_TYPE>
using OpMixScalarTimesDiv = FormsIntegrators<DomainEleOp>::Assembly<
A>::BiLinearForm<I>::OpMixScalarTimesDiv<SPACE_DIM, COORD_TYPE>;
// Flux is applied by Lagrange Multiplie
using BoundaryNaturalBC = NaturalBC<BoundaryEleOp>::Assembly<A>::LinearForm<I>;
using OpFluidFlux =
BoundaryNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 1>;
FTensor::Index<'i', SPACE_DIM> i;
FTensor::Index<'j', SPACE_DIM> j;
FTensor::Index<'k', SPACE_DIM> k;
FTensor::Index<'l', SPACE_DIM> l;
constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
// 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.03; // 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) {
if (coord_type == CYLINDRICAL)
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) {
FTensor::Ddg<double, SPACE_DIM, SPACE_DIM> t_D;
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) {
MoFEMFunctionBegin;
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);
}
MoFEMFunctionReturn(0);
};
/**
* @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());
Range r;
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());
Range r;
r.insert_list(ents_vec.begin(), it);
return r;
};
#include <FreeSurfaceOps.hpp>
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
*/
MoFEMErrorCode tsSetUp(TS ts);
SmartPetscObj<VecScatter> getScatter(Vec x, Vec y, enum FR fr);
SmartPetscObj<Vec> getSubVector();
SmartPetscObj<DM> solverSubDM;
SmartPetscObj<Vec> globSol;
FreeSurface *fsRawPtr;
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) {}
MoFEMErrorCode runProblem();
MoFEMErrorCode makeRefProblem();
MoFEM::Interface &mField;
private:
MoFEMErrorCode readMesh();
MoFEMErrorCode setupProblem();
MoFEMErrorCode boundaryCondition();
MoFEMErrorCode projectData();
MoFEMErrorCode assembleSystem();
MoFEMErrorCode solveSystem();
/// @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
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);
friend struct TSPrePostProc;
};
//! [Run programme]
MoFEMErrorCode FreeSurface::runProblem() {
MoFEMFunctionBegin;
CHKERR readMesh();
CHKERR setupProblem();
CHKERR boundaryCondition();
CHKERR assembleSystem();
CHKERR solveSystem();
MoFEMFunctionReturn(0);
}
//! [Run programme]
//! [Read mesh]
MoFEMErrorCode FreeSurface::readMesh() {
MoFEMFunctionBegin;
MOFEM_LOG("FS", Sev::inform) << "Read mesh for problem";
auto simple = mField.getInterface<Simple>();
simple->getParentAdjacencies() = true;
simple->getBitRefLevel() = BitRefLevel();
CHKERR simple->getOptions();
CHKERR simple->loadFile();
MoFEMFunctionReturn(0);
}
//! [Read mesh]
//! [Set up problem]
MoFEMErrorCode FreeSurface::setupProblem() {
MoFEMFunctionBegin;
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,
LAST_COORDINATE_SYSTEM, &coord_type, PETSC_NULL);
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 simple = mField.getInterface<Simple>();
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);
CHKERR PetscOptionsGetScalar(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
CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, U_FIELD_DIM);
// 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,
U_FIELD_DIM);
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 = [&]() {
MoFEMFunctionBegin;
// 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();
MoFEMFunctionReturn(0);
};
CHKERR set_problem_bit();
CHKERR simple->setUp();
MoFEMFunctionReturn(0);
}
//! [Set up problem]
//! [Boundary condition]
MoFEMErrorCode FreeSurface::boundaryCondition() {
MoFEMFunctionBegin;
#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 simple = mField.getInterface<Simple>();
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 = [&]() {
MoFEMFunctionBegin;
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);
MoFEMFunctionReturn(0);
};
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(
new DomianParentEle(mField));
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) {
MoFEMFunctionBegin;
auto &pip = fe->getOpPtrVector();
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
CHKERR setParentDofs(
fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
new DomianParentEle(mField));
AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
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(
new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
CHKERR add_parent_field(fe, UDO::OPROW, "G");
pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
pip.push_back(
new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
MoFEMFunctionReturn(0);
};
auto post_proc = [&](auto exe_test) {
MoFEMFunctionBegin;
auto post_proc_fe = boost::make_shared<PostProcEleDomain>(mField);
post_proc_fe->exeTestHook = exe_test;
CHKERR set_generic(post_proc_fe);
using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
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");
MoFEMFunctionReturn(0);
};
auto solve_init = [&](auto exe_test) {
MoFEMFunctionBegin;
pip_mng->getOpDomainRhsPipeline().clear();
pip_mng->getOpDomainLhsPipeline().clear();
auto set_domain_rhs = [&](auto fe) {
MoFEMFunctionBegin;
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));
MoFEMFunctionReturn(0);
};
auto set_domain_lhs = [&](auto fe) {
MoFEMFunctionBegin;
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));
MoFEMFunctionReturn(0);
};
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) {
MoFEMFunctionBegin;
CHKERR SNESMonitorSet(solver,
(MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
void *))MoFEMSNESMonitorFields,
(void *)(snes_ctx_ptr.get()), nullptr);
MoFEMFunctionReturn(0);
};
auto print_section_field = [&]() {
MoFEMFunctionBegin;
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);
}
MoFEMFunctionReturn(0);
};
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);
MoFEMFunctionReturn(0);
};
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);
MoFEMFunctionReturn(0);
}
//! [Boundary condition]
//! [Data projection]
MoFEMErrorCode FreeSurface::projectData() {
MoFEMFunctionBegin;
// FIXME: Functionality in this method should be improved (projection should
// be done field by field), generalised, and move to become core
// functionality.
auto simple = mField.getInterface<Simple>();
auto pip_mng = mField.getInterface<PipelineManager>();
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(
new DomianParentEle(mField));
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(
new BoundaryParentEle(mField));
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) {
MoFEMFunctionBegin;
auto set_domain_rhs = [&](auto fe) {
MoFEMFunctionBegin;
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();
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
CHKERR setParentDofs(
eval_fe_ptr, "", UDO::OPSPACE, bit(get_skin_projection_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
new DomianParentEle(mField));
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",
bit(get_skin_projection_bit()));
pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW, "P",
bit(get_skin_projection_bit()));
pip.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW, "H",
bit(get_skin_projection_bit()));
pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW, "G",
bit(get_skin_projection_bit()));
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,
bit(get_projection_bit())));
auto assemble_fe_ptr = boost::make_shared<DomianParentEle>(mField);
{
auto &pip = assemble_fe_ptr->getOpPtrVector();
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
CHKERR setParentDofs(
assemble_fe_ptr, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
new DomianParentEle(mField));
return fe_parent;
},
QUIET, Sev::noisy);
CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "U",
bit(get_skin_parent_bit()));
pip.push_back(new OpDomainAssembleVector("U", u_ptr));
CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "P",
bit(get_skin_parent_bit()));
pip.push_back(new OpDomainAssembleScalar("P", p_ptr));
CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "H",
bit(get_skin_parent_bit()));
pip.push_back(new OpDomainAssembleScalar("H", h_ptr));
CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "G",
bit(get_skin_parent_bit()));
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())));
MoFEMFunctionReturn(0);
};
auto set_domain_lhs = [&](auto fe) {
MoFEMFunctionBegin;
auto &pip = fe->getOpPtrVector();
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
CHKERR setParentDofs(
fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
new DomianParentEle(mField));
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",
bit(get_skin_parent_bit()));
CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U",
bit(get_skin_parent_bit()));
pip.push_back(new OpDomainMassU("U", "U"));
CHKERR add_parent_field_domain(fe, UDO::OPROW, "P",
bit(get_skin_parent_bit()));
CHKERR add_parent_field_domain(fe, UDO::OPCOL, "P",
bit(get_skin_parent_bit()));
pip.push_back(new OpDomainMassP("P", "P"));
CHKERR add_parent_field_domain(fe, UDO::OPROW, "H",
bit(get_skin_parent_bit()));
CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H",
bit(get_skin_parent_bit()));
pip.push_back(new OpDomainMassH("H", "H"));
CHKERR add_parent_field_domain(fe, UDO::OPROW, "G",
bit(get_skin_parent_bit()));
CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G",
bit(get_skin_parent_bit()));
pip.push_back(new OpDomainMassG("G", "G"));
MoFEMFunctionReturn(0);
};
auto set_bdy_rhs = [&](auto fe) {
MoFEMFunctionBegin;
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();
CHKERR setParentDofs(
eval_fe_ptr, "", UDO::OPSPACE, bit(get_skin_projection_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
new BoundaryParentEle(mField));
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",
bit(get_skin_projection_bit()));
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,
bit(get_projection_bit())));
auto assemble_fe_ptr = boost::make_shared<BoundaryParentEle>(mField);
{
auto &pip = assemble_fe_ptr->getOpPtrVector();
CHKERR setParentDofs(
assemble_fe_ptr, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
new BoundaryParentEle(mField));
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 &) {
MoFEMFunctionBegin;
if (lPtr->size() != getGaussPts().size2()) {
lPtr->resize(getGaussPts().size2());
lPtr->clear();
}
MoFEMFunctionReturn(0);
}
private:
boost::shared_ptr<VectorDouble> lPtr;
};
pip.push_back(new OpLSize(l_ptr));
CHKERR add_parent_field_bdy(assemble_fe_ptr, UDO::OPROW, "L",
bit(get_skin_parent_bit()));
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())));
MoFEMFunctionReturn(0);
};
auto set_bdy_lhs = [&](auto fe) {
MoFEMFunctionBegin;
auto &pip = fe->getOpPtrVector();
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);
// 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",
bit(get_skin_parent_bit()));
CHKERR add_parent_field_bdy(fe, UDO::OPCOL, "L",
bit(get_skin_parent_bit()));
pip.push_back(new OpBoundaryMassL("L", "L"));
MoFEMFunctionReturn(0);
};
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>();
MoFEMFunctionBegin;
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,
mField.get_comm_size());
// set ghost nodes
CHKERR prb_mng->partitionGhostDofsOnDistributedMesh("SUB_SOLVER");
MoFEMFunctionReturn(0);
};
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";
return SmartPetscObj<DM>();
};
auto create_dummy_dm = [&]() {
auto dummy_dm = createDM(mField.get_comm(), "DMMOFEM");
CHK_THROW_MESSAGE(DMMoFEMCreateMoFEM(dummy_dm, &mField,
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) {
MoFEMFunctionBegin;
for (auto &v : ent_ptr->getEntFieldData()) {
v = 0;
}
MoFEMFunctionReturn(0);
};
auto solve = [&](auto S) {
MoFEMFunctionBegin;
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);
MoFEMFunctionReturn(0);
};
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 = [&]() {
MoFEMFunctionBegin;
int n, ystart;
CHKERR VecGetLocalSize(v, &n);
CHKERR VecGetOwnershipRange(v, &ystart, NULL);
CHKERR ISCreateStride(PETSC_COMM_SELF, n, ystart, 1, &iy);
MoFEMFunctionReturn(0);
};
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);
CHK_THROW_MESSAGE(
DMSubDomainRestrict(simple->getDM(), s, PETSC_NULL, dummy_dm),
"restrict");
Vec X0, Xdot;
CHK_THROW_MESSAGE(DMGetNamedGlobalVector(dummy_dm, "TSTheta_X0", &X0),
"get X0");
CHK_THROW_MESSAGE(
DMGetNamedGlobalVector(dummy_dm, "TSTheta_Xdot", &Xdot),
"get Xdot");
auto forward_ghost = [](auto g) {
MoFEMFunctionBegin;
CHKERR VecGhostUpdateBegin(g, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(g, INSERT_VALUES, SCATTER_FORWARD);
MoFEMFunctionReturn(0);
};
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()) {
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);
}
MoFEMFunctionReturn(0);
};
// postprocessing (only for debugging proposes)
auto post_proc = [&](auto exe_test) {
MoFEMFunctionBegin;
auto post_proc_fe = boost::make_shared<PostProcEleDomain>(mField);
auto &pip = post_proc_fe->getOpPtrVector();
post_proc_fe->exeTestHook = exe_test;
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
CHKERR setParentDofs(
post_proc_fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
new DomianParentEle(mField));
AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
fe_parent->getOpPtrVector(), {H1});
return fe_parent;
},
QUIET, Sev::noisy);
using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
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",
bit(get_skin_parent_bit()));
pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "P",
bit(get_skin_parent_bit()));
pip.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "H",
bit(get_skin_parent_bit()));
pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "G",
bit(get_skin_parent_bit()));
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");
MoFEMFunctionReturn(0);
};
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());
MoFEMFunctionReturn(0);
}
//! [Data projection]
//! [Push operators to pip]
MoFEMErrorCode FreeSurface::assembleSystem() {
MoFEMFunctionBegin;
auto simple = mField.getInterface<Simple>();
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(
new DomianParentEle(mField));
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(
new BoundaryParentEle(mField));
return fe_parent;
},
QUIET, Sev::noisy);
};
auto test_bit_child = [](FEMethod *fe_ptr) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
get_start_bit() + nb_levels - 1);
};
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) {
MoFEMFunctionBegin;
auto &pip = fe->getOpPtrVector();
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
CHKERR setParentDofs(
fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
new DomianParentEle(mField));
AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
fe_parent->getOpPtrVector(), {H1});
return fe_parent;
},
QUIET, Sev::noisy);
CHKERR add_parent_field_domain(fe, UDO::OPROW, "U");
pip.push_back(
new OpCalculateVectorFieldValuesDot<U_FIELD_DIM>("U", dot_u_ptr));
pip.push_back(new OpCalculateVectorFieldValues<U_FIELD_DIM>("U", u_ptr));
pip.push_back(new OpCalculateVectorFieldGradient<U_FIELD_DIM, SPACE_DIM>(
"U", grad_u_ptr));
switch (coord_type) {
case CARTESIAN:
pip.push_back(
new OpCalculateDivergenceVectorFieldValues<SPACE_DIM, CARTESIAN>(
"U", div_u_ptr));
break;
case CYLINDRICAL:
pip.push_back(
new OpCalculateDivergenceVectorFieldValues<SPACE_DIM, CYLINDRICAL>(
"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(
new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
CHKERR add_parent_field_domain(fe, UDO::OPROW, "G");
pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
pip.push_back(
new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
CHKERR add_parent_field_domain(fe, UDO::OPROW, "P");
pip.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
MoFEMFunctionReturn(0);
};
auto set_domain_rhs = [&](auto fe) {
MoFEMFunctionBegin;
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);
}));
MoFEMFunctionReturn(0);
};
auto set_domain_lhs = [&](auto fe) {
MoFEMFunctionBegin;
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;
case CYLINDRICAL:
pip.push_back(new OpMixScalarTimesDiv<CYLINDRICAL>(
"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);
}));
}
MoFEMFunctionReturn(0);
};
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) {
MoFEMFunctionBegin;
auto &pip = fe->getOpPtrVector();
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);
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;
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
op_bdy_side->getOpPtrVector(), {H1});
CHKERR setParentDofs(
op_bdy_side->getSideFEPtr(), "", UDO::OPSPACE,
bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
new DomianParentEle(mField));
AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
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(
new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
// 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()));
}
}
MoFEMFunctionReturn(0);
};
auto set_boundary_lhs = [&](auto fe) {
MoFEMFunctionBegin;
auto &pip = fe->getOpPtrVector();
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);
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;
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
op_bdy_side->getOpPtrVector(), {H1});
CHKERR setParentDofs(
op_bdy_side->getSideFEPtr(), "", UDO::OPSPACE,
bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
new DomianParentEle(mField));
AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
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(
new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
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()));
}
}
MoFEMFunctionReturn(0);
};
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;
MoFEMFunctionReturn(0);
}
//! [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 {
Monitor(
SmartPetscObj<DM> dm, boost::shared_ptr<moab::Core> post_proc_mesh,
boost::shared_ptr<PostProcEleDomainCont> post_proc,
boost::shared_ptr<PostProcEleBdyCont> post_proc_edge,
std::pair<boost::shared_ptr<BoundaryEle>, boost::shared_ptr<VectorDouble>>
p)
: dM(dm), postProcMesh(post_proc_mesh), postProc(post_proc),
postProcEdge(post_proc_edge), liftFE(p.first), liftVec(p.second) {}
MoFEMErrorCode postProcess() {
MoFEMFunctionBegin;
MOFEM_LOG("FS", Sev::verbose) << "Monitor";
constexpr int save_every_nth_step = 1;
if (ts_step % save_every_nth_step == 0) {
MOFEM_LOG("FS", Sev::verbose) << "Mesh pre proc";
MoFEM::Interface *m_field_ptr;
CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
auto post_proc_begin =
boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
postProcMesh);
auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
*m_field_ptr, postProcMesh);
CHKERR DMoFEMPreProcessFiniteElements(dM, post_proc_begin->getFEMethod());
CHKERR DMoFEMLoopFiniteElements(dM, "dFE", postProc,
this->getCacheWeakPtr());
CHKERR DMoFEMLoopFiniteElements(dM, "bFE", postProcEdge,
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();
CHKERR DMoFEMLoopFiniteElements(dM, "bFE", liftFE, this->getCacheWeakPtr());
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];
MoFEMFunctionReturn(0);
}
private:
SmartPetscObj<DM> dM;
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]
MoFEMErrorCode FreeSurface::solveSystem() {
MoFEMFunctionBegin;
using UDO = ForcesAndSourcesCore::UserDataOperator;
auto simple = mField.getInterface<Simple>();
auto pip_mng = mField.getInterface<PipelineManager>();
auto create_solver_dm = [&](auto dm) -> SmartPetscObj<DM> {
DM subdm;
auto setup_subdm = [&](auto dm) {
MoFEMFunctionBegin;
auto simple = mField.getInterface<Simple>();
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);
MoFEMFunctionReturn(0);
};
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(
new DomianParentEle(mField));
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(
get_start_bit() + nb_levels - 1);
};
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>();
AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
post_proc_fe->getOpPtrVector(), {H1});
CHKERR setParentDofs(
post_proc_fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
[&]() {
boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
new DomianParentEle(mField));
AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
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(
new OpCalculateVectorFieldValues<U_FIELD_DIM>("U", u_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateVectorFieldGradient<U_FIELD_DIM, SPACE_DIM>("U",
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(
new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
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(
new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
using OpPPMap = OpPostProcMapInMoab<2, 2>;
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(
get_start_bit() + nb_levels - 1);
};
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) {}
MoFEMErrorCode doWork(int side, EntityType type,
EntitiesFieldData::EntData &data) {
MoFEMFunctionBegin;
auto t_l = getFTensor0FromVec(*ptrL);
auto t_n_fe = getFTensor1NormalsAtGaussPts();
ptrNormal->resize(SPACE_DIM, getGaussPts().size2(), false);
auto t_n = getFTensor1FromMat<SPACE_DIM>(*ptrNormal);
for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
t_n(i) = t_n_fe(i) * t_l / std::sqrt(t_n_fe(i) * t_n_fe(i));
++t_n_fe;
++t_l;
++t_n;
}
MoFEMFunctionReturn(0);
};
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(
new OpCalculateVectorFieldValues<U_FIELD_DIM>("U", u_ptr));
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));
auto op_ptr = new BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE);
op_ptr->doWorkRhsHook = [&](DataOperator *base_op_ptr, int side,
EntityType type,
EntitiesFieldData::EntData &data) {
MoFEMFunctionBegin;
MoFEMFunctionReturn(0);
};
using OpPPMap = OpPostProcMapInMoab<2, 2>;
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(
get_start_bit() + nb_levels - 1);
};
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) {
MoFEMFunctionBegin;
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);
MoFEMFunctionReturn(0);
};
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
CHKERR TSSetFromOptions(ts);
CHKERR ptr->tsSetUp(ts);
CHKERR TSSetUp(ts);
auto print_fields_in_section = [&]() {
MoFEMFunctionBegin;
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);
}
MoFEMFunctionReturn(0);
};
CHKERR print_fields_in_section();
CHKERR TSSolve(ts, ptr->globSol);
}
MoFEMFunctionReturn(0);
}
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";
MoFEM::Core::Initialize(&argc, &argv, param_file, help);
// 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";
CHKERR DMRegister_MoFEM(dm_name);
//! [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]
}
CATCH_ERRORS;
CHKERR MoFEM::Core::Finalize();
#ifdef PYTHON_INIT_SURFACE
if (Py_FinalizeEx() < 0) {
exit(120);
}
#endif
}
std::vector<Range>
FreeSurface::findEntitiesCrossedByPhaseInterface(size_t overlap) {
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_LOG_CHANNEL("SYNC");
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;
MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::noisy);
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) {
MoFEMFunctionBegin;
CHKERR bit_mng->updateRangeByChildren(p, c);
c = c.subset_by_dimension(SPACE_DIM);
MoFEMFunctionReturn(0);
};
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) {
Range l;
CHK_THROW_MESSAGE(
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;
}
MoFEMErrorCode FreeSurface::refineMesh(size_t overlap) {
MoFEMFunctionBegin;
auto bit_mng = mField.getInterface<BitRefManager>();
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*/) {
MoFEMFunctionBegin;
// 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);
}
MoFEMFunctionReturn(0);
};
// resolve skin between refined levels
auto set_skins = [&]() {
MoFEMFunctionBegin;
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;
CHK_THROW_MESSAGE(
mField.getInterface<BitRefManager>()->getEntitiesByDimAndRefLevel(
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);
MoFEMFunctionReturn(0);
};
// take last level, remove childs on boarder, and set bit
auto set_current = [&]() {
MoFEMFunctionBegin;
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);
MoFEMFunctionReturn(0);
};
// set bits to levels
CHKERR set_levels(findEntitiesCrossedByPhaseInterface(overlap));
// 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");
}
}
MoFEMFunctionReturn(0);
};
MoFEMErrorCode FreeSurface::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) {
MoFEMFunctionBegin;
/**
* @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(
new OpAddParentEntData(
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(
new OpAddParentEntData(
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),
nb_levels);
MoFEMFunctionReturn(0);
}
MoFEMErrorCode TSPrePostProc::tsPreProc(TS ts) {
MoFEMFunctionBegin;
if (auto ptr = tsPrePostProc.lock()) {
/**
* @brief cut-off values at nodes, i.e. abs("H") <= 1
*
*/
auto cut_off_dofs = [&]() {
MoFEMFunctionBegin;
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) {
MoFEMFunctionBegin;
for (auto &h : ent_ptr->getEntFieldData()) {
h = cut_off(h);
}
MoFEMFunctionReturn(0);
};
auto field_blas = m_field.getInterface<FieldBlas>();
CHKERR field_blas->fieldLambdaOnEntities(cut_off_verts, "H",
&current_verts);
MoFEMFunctionReturn(0);
};
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>();
// 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 = [&]() {
MoFEMFunctionBegin;
MOFEM_LOG("FS", Sev::inform) << "Refine problem";
CHKERR ptr->fsRawPtr->refineMesh(refine_overlap);
CHKERR ptr->fsRawPtr->projectData();
MoFEMFunctionReturn(0);
};
// set new jacobin operator, since problem and thus tangent matrix size has
// changed
auto set_jacobian_operators = [&]() {
MoFEMFunctionBegin;
ptr->subB = createDMMatrix(ptr->solverSubDM);
CHKERR KSPReset(ptr->subKSP);
MoFEMFunctionReturn(0);
};
// set new solution
auto set_solution = [&]() {
MoFEMFunctionBegin;
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);
MoFEMFunctionReturn(0);
};
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_LOG_CHANNEL("SYNC");
MOFEM_TAG_AND_LOG("SYNC", Sev::verbose, "FS") << "PreProc done";
MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::verbose);
}
MoFEMFunctionReturn(0);
}
MoFEMErrorCode TSPrePostProc::tsPostProc(TS ts) {
if (auto ptr = tsPrePostProc.lock()) {
auto &m_field = ptr->fsRawPtr->mField;
MOFEM_LOG_CHANNEL("SYNC");
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) {
MoFEMFunctionBegin;
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) {
MoFEMFunctionBegin;
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);
MoFEMFunctionReturn(0);
};
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);
}
MoFEMFunctionReturn(0);
}
MoFEMErrorCode TSPrePostProc::tsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t,
PetscReal a, Mat A, Mat B,
void *ctx) {
MoFEMFunctionBegin;
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) {
MoFEMFunctionBegin;
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);
MoFEMFunctionReturn(0);
};
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());
}
MoFEMFunctionReturn(0);
}
MoFEMErrorCode TSPrePostProc::tsMonitor(TS ts, PetscInt step, PetscReal t,
Vec u, void *ctx) {
MoFEMFunctionBegin;
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());
}
MoFEMFunctionReturn(0);
}
MoFEMErrorCode TSPrePostProc::pcSetup(PC pc) {
MoFEMFunctionBegin;
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);
}
MoFEMFunctionReturn(0);
};
MoFEMErrorCode TSPrePostProc::pcApply(PC pc, Vec pc_f, Vec pc_x) {
MoFEMFunctionBegin;
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);
}
MoFEMFunctionReturn(0);
};
SmartPetscObj<VecScatter> TSPrePostProc::getScatter(Vec x, Vec y, enum FR fr) {
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");
}
return SmartPetscObj<VecScatter>(s);
}
}
CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "No prb pinter");
return SmartPetscObj<VecScatter>();
}
SmartPetscObj<Vec> TSPrePostProc::getSubVector() {
return createDMVector(solverSubDM);
}
MoFEMErrorCode TSPrePostProc::tsSetUp(TS ts) {
auto &m_field = fsRawPtr->mField;
auto simple = m_field.getInterface<Simple>();
MoFEMFunctionBegin;
auto dm = simple->getDM();
globRes = vectorDuplicate(globSol);
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) {
MoFEMFunctionBegin;
CHKERR SNESMonitorSet(snes,
(MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
void *))MoFEMSNESMonitorFields,
(void *)(snes_ctx_ptr.get()), nullptr);
MoFEMFunctionReturn(0);
};
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);
snesCtxPtr = getDMSnesCtx(solverSubDM);
tsCtxPtr = getDMTsCtx(solverSubDM);
CHKERR TSSetPreStep(ts, tsPreProc);
CHKERR TSSetPostStep(ts, tsPostProc);
MoFEMFunctionReturn(0);
}
NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEM::NaturalBC::Assembly::LinearForm
Definition: Natural.hpp:67
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:589
get_D
auto get_D
Definition: free_surface.cpp:254
rho_ave
double rho_ave
Definition: free_surface.cpp:180
TSPrePostProc::solverSubDM
SmartPetscObj< DM > solverSubDM
Definition: free_surface.cpp:421
MoFEM::OpPostProcMapInMoab::DataMapMat
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
Definition: PostProcBrokenMeshInMoabBase.hpp:712
FreeSurface::findEntitiesCrossedByPhaseInterface
std::vector< Range > findEntitiesCrossedByPhaseInterface(size_t overlap)
Find entities on refinement levels.
Definition: free_surface.cpp:2453
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
get_fe_bit
auto get_fe_bit
Definition: free_surface.cpp:316
g
constexpr double g
Definition: shallow_wave.cpp:63
FreeSurface::projectData
MoFEMErrorCode projectData()
[Boundary condition]
Definition: free_surface.cpp:953
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
MoFEM::debug
const static bool debug
Definition: ConstrainMatrixCtx.cpp:9
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
FreeSurfaceOps::OpNormalForceRhs
Definition: FreeSurfaceOps.hpp:92
MoFEM::DMMoFEMAddSubFieldCol
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:280
MoFEM::DataOperator
base operator to do operations at Gauss Pt. level
Definition: DataOperators.hpp:24
MOFEM_LOG_SEVERITY_SYNC
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
Definition: LogManager.hpp:352
kernel_oscillation
auto kernel_oscillation
Definition: free_surface.cpp:262
Monitor::liftFE
boost::shared_ptr< BoundaryEle > liftFE
Definition: free_surface.cpp:2039
CARTESIAN
@ CARTESIAN
Definition: definitions.h:128
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
EntityHandle
get_M_dh
auto get_M_dh
Definition: free_surface.cpp:252
MoFEM::OpCalculateDivergenceVectorFieldValues
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:488
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
FreeSurfaceOps::OpWettingAngleRhs
Definition: FreeSurfaceOps.hpp:138
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:373
get_skin_child_bit
auto get_skin_child_bit
Definition: free_surface.cpp:155
MoFEM::OpFluxRhsImpl
Definition: Natural.hpp:39
cut_off
auto cut_off
Definition: free_surface.cpp:205
get_M
auto get_M
Definition: free_surface.cpp:251
MoFEM::OpRunParent
Operator to execute finite element instance on parent element. This operator is typically used to pro...
Definition: MeshProjectionDataOperators.hpp:17
MoFEM::PipelineManager::ElementsAndOpsByDim
Definition: PipelineManager.hpp:38
MoFEM::createTS
auto createTS(MPI_Comm comm)
Definition: PetscSmartObj.hpp:249
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::TsSetIFunction
PetscErrorCode TsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t, Vec F, void *ctx)
Set IFunction for TS solver.
Definition: TsCtx.cpp:56
get_skin_parent_bit
auto get_skin_parent_bit
Definition: free_surface.cpp:154
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:456
init_h
auto init_h
Initialisation function.
Definition: free_surface.cpp:297
get_f_dh
auto get_f_dh
Definition: free_surface.cpp:224
mu_m
double mu_m
Definition: free_surface.cpp:164
help
static char help[]
Definition: activate_deactivate_dofs.cpp:13
FreeSurfaceOps::OpWettingAngleLhs
Definition: FreeSurfaceOps.hpp:246
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
FreeSurface::assembleSystem
MoFEMErrorCode assembleSystem()
[Data projection]
Definition: free_surface.cpp:1581
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:468
MoFEM::getPetscObject
PetscObject getPetscObject(T obj)
Definition: PetscSmartObj.hpp:9
d_phase_function_h
auto d_phase_function_h
Definition: free_surface.cpp:217
my_min
auto my_min
Definition: free_surface.cpp:204
nb_levels
constexpr int nb_levels
Definition: level_set.cpp:58
MoFEM::CoreInterface::get_field_bit_number
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
PostProcEleDomain
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
Definition: dynamic_first_order_con_law.cpp:52
md
double md
Definition: free_surface.cpp:176
MoFEM::OpCalculateScalarFieldValuesDot
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
Definition: UserDataOperators.hpp:273
MoFEM::ForcesAndSourcesCore::UserDataOperator::OpType
OpType
Controls loop over entities on element.
Definition: ForcesAndSourcesCore.hpp:566
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
FreeSurface::runProblem
MoFEMErrorCode runProblem()
[Run programme]
Definition: free_surface.cpp:524
MoFEM.hpp
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
MoFEM::DMoFEMMeshToLocalVector
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:523
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
get_f
auto get_f
Definition: free_surface.cpp:223
MoFEM::createKSP
auto createKSP(MPI_Comm comm)
Definition: PetscSmartObj.hpp:261
OpBoundaryAssembleScalar
FormsIntegrators< BoundaryEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpBoundaryAssembleScalar
Definition: free_surface.cpp:123
W
double W
Definition: free_surface.cpp:168
FreeSurfaceOps::OpLhsH_dG
Lhs for H dH.
Definition: FreeSurfaceOps.hpp:1117
OpBase
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
tsPrePostProc
static boost::weak_ptr< TSPrePostProc > tsPrePostProc
Definition: dynamic_first_order_con_law.cpp:404
BASE_DIM
constexpr int BASE_DIM
Definition: dg_projection.cpp:15
test_bit_child
auto test_bit_child
lambda function used to select elements on which finite element pipelines are executed.
Definition: hanging_node_approx.cpp:173
Monitor::postProcess
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: dynamic_first_order_con_law.cpp:789
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
TSPrePostProc::getScatter
SmartPetscObj< VecScatter > getScatter(Vec x, Vec y, enum FR fr)
Definition: free_surface.cpp:3107
get_M3
auto get_M3
Definition: free_surface.cpp:235
sdf.r
int r
Definition: sdf.py:8
MoFEM::createDMMatrix
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:1056
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
order
constexpr int order
Definition: dg_projection.cpp:18
capillary_tube
auto capillary_tube
Definition: free_surface.cpp:280
OpMixScalarTimesDiv
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMixScalarTimesDiv< SPACE_DIM, COORD_TYPE > OpMixScalarTimesDiv
Definition: free_surface.cpp:127
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1319
I
constexpr IntegrationType I
Definition: operators_tests.cpp:31
MoFEM::CoreInterface::get_dofs
virtual const DofEntity_multiIndex * get_dofs() const =0
Get the dofs object.
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::TsMonitorSet
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition: TsCtx.cpp:263
FreeSurfaceOps::OpRhsH
Definition: FreeSurfaceOps.hpp:795
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl
Approximate field values for given petsc vector.
Definition: UserDataOperators.hpp:597
cylindrical
auto cylindrical
Definition: free_surface.cpp:189
OpDomainAssembleVector
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesVector< BASE_DIM, SPACE_DIM, 1 > OpDomainAssembleVector
Definition: free_surface.cpp:119
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
MoFEM::DMMoFEMSetDestroyProblem
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition: DMMoFEM.cpp:438
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:677
U_FIELD_DIM
constexpr int U_FIELD_DIM
Definition: free_surface.cpp:81
OpBoundaryMassL
FormsIntegrators< BoundaryEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, 1 > OpBoundaryMassL
Definition: free_surface.cpp:116
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:178
TSPrePostProc
Set of functions called by PETSc solver used to refine and update mesh.
Definition: dynamic_first_order_con_law.cpp:383
save_every_nth_step
int save_every_nth_step
Definition: photon_diffusion.cpp:67
kernel_eye
auto kernel_eye
Definition: free_surface.cpp:272
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
FreeSurfaceOps::OpLhsH_dU
Lhs for H dU.
Definition: FreeSurfaceOps.hpp:915
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
TSPrePostProc::tsMonitor
static MoFEMErrorCode tsMonitor(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Wrapper for TS monitor.
Definition: free_surface.cpp:3052
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
wetting_angle
auto wetting_angle
Definition: free_surface.cpp:312
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
FreeSurfaceOps::OpLoopSideGetDataForSideEle
Definition: FreeSurfaceOps.hpp:576
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
MoFEM::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition: PetscSmartObj.hpp:141
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition: seepage.cpp:57
BoundaryEleOp
R
@ R
Definition: free_surface.cpp:396
MoFEM::BitRefManager::getEntitiesByDimAndRefLevel
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
Definition: BitRefManager.cpp:822
rho_diff
double rho_diff
Definition: free_surface.cpp:181
MoFEM::get_temp_meshset_ptr
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition: Templates.hpp:1886
get_current_bit
auto get_current_bit
dofs bit used to do calculations
Definition: free_surface.cpp:151
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
MoFEM::DMMoFEMCreateSubDM
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:215
double
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
convert.type
type
Definition: convert.py:64
OpDomainMassH
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, 1 > OpDomainMassH
Definition: free_surface.cpp:112
get_global_size
auto get_global_size
Definition: free_surface.cpp:320
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:317
FreeSurfaceOps::OpLhsH_dH
Lhs for H dH.
Definition: FreeSurfaceOps.hpp:973
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
eta2
double eta2
Definition: free_surface.cpp:173
TSPrePostProc::tsSetIFunction
static MoFEMErrorCode tsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t, Vec f, void *ctx)
Definition: free_surface.cpp:2997
TSPrePostProc::tsSetUp
MoFEMErrorCode tsSetUp(TS ts)
Used to setup TS solver.
Definition: free_surface.cpp:3131
MoFEM::Problem::getSubData
boost::shared_ptr< SubProblemData > & getSubData() const
Get main problem of sub-problem is.
Definition: ProblemsMultiIndices.hpp:119
MoFEM::PostProcBrokenMeshInMoabBaseContImpl
Definition: PostProcBrokenMeshInMoabBase.hpp:890
PlasticOps::Monitor
Definition: PlasticOpsMonitor.hpp:9
h
double h
Definition: photon_diffusion.cpp:60
MoFEM::FieldEntity::getLoBitNumberUId
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
Definition: FieldEntsMultiIndices.hpp:222
OpDomainMassP
OpDomainMassH OpDomainMassP
Definition: free_surface.cpp:113
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
get_projection_bit
auto get_projection_bit
Definition: free_surface.cpp:156
BoundaryEleOp
BoundaryEle::UserDataOperator BoundaryEleOp
Definition: child_and_parent.cpp:40
get_M0_dh
auto get_M0_dh
Definition: free_surface.cpp:227
TSPrePostProc::tsPostProc
static MoFEMErrorCode tsPostProc(TS ts)
Post process time step.
Definition: free_surface.cpp:2987
MoFEM::createPC
auto createPC(MPI_Comm comm)
Definition: PetscSmartObj.hpp:267
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
SideEle
ElementsAndOps< SPACE_DIM >::SideEle SideEle
Definition: plastic.cpp:61
eta
double eta
Definition: free_surface.cpp:172
MoFEM::CoreInterface::get_comm_size
virtual int get_comm_size() const =0
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
SideOp
SideEle::UserDataOperator SideOp
Definition: free_surface.cpp:97
FreeSurfaceOps::OpRhsG
Definition: FreeSurfaceOps.hpp:1181
phase_function
auto phase_function
Definition: free_surface.cpp:213
wetting_angle_sub_stepping
auto wetting_angle_sub_stepping
Definition: free_surface.cpp:196
MoFEM::solve
static PetscErrorCode solve(Mat mat, Vec x, Vec y)
Definition: Schur.cpp:1296
Monitor::postProcMesh
boost::shared_ptr< moab::Core > postProcMesh
Definition: free_surface.cpp:2036
MoFEM::TsSetIJacobian
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:169
Monitor::dM
SmartPetscObj< DM > dM
Definition: dynamic_first_order_con_law.cpp:873
MoFEM::DMoFEMPreProcessFiniteElements
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:546
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
MoFEM::LogManager::SeverityLevel
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
MoFEM::DMoFEMPostProcessFiniteElements
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:556
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
MoFEM::AddFluxToRhsPipelineImpl
Definition: Natural.hpp:46
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
debug
static const bool debug
Definition: dm_create_subdm.cpp:12
MoFEM::DMMoFEMCreateMoFEM
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:114
MoFEM::ForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: ForcesAndSourcesCore.hpp:482
TSPrePostProc::snesCtxPtr
boost::shared_ptr< SnesCtx > snesCtxPtr
Definition: free_surface.cpp:464
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::DataOperator::doWorkRhsHook
DoWorkRhsHookFunType doWorkRhsHook
Definition: DataOperators.hpp:62
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:139
FreeSurfaceOps::OpLhsU_dG
Lhs for G dH.
Definition: FreeSurfaceOps.hpp:732
COORD_TYPE
constexpr CoordinateTypes COORD_TYPE
Definition: tensor_divergence_operator.cpp:36
BiLinearForm
get_skin_projection_bit
auto get_skin_projection_bit
Definition: free_surface.cpp:157
main
int main(int argc, char *argv[])
Definition: activate_deactivate_dofs.cpp:15
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', SPACE_DIM >
MoFEM::IntegrationType
IntegrationType
Form integrator integration types.
Definition: FormsIntegrators.hpp:136
mu_diff
double mu_diff
Definition: free_surface.cpp:183
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:417
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1561
convert.n
n
Definition: convert.py:82
rho_p
double rho_p
Definition: free_surface.cpp:165
FreeSurfaceOps::OpLhsG_dG
Lhs for H dH.
Definition: FreeSurfaceOps.hpp:1321
get_M3_dh
auto get_M3_dh
Definition: free_surface.cpp:244
MoFEM::Problem::numeredRowDofsPtr
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredRowDofsPtr
store DOFs on rows for this problem
Definition: ProblemsMultiIndices.hpp:73
get_M2_dh
auto get_M2_dh
Definition: free_surface.cpp:233
kappa
double kappa
Definition: free_surface.cpp:185
MoFEM::TSMethod::ts_step
PetscInt ts_step
time step number
Definition: LoopMethods.hpp:163
FreeSurfaceOps.hpp
rho_m
double rho_m
Definition: free_surface.cpp:163
FreeSurface
Definition: free_surface.cpp:471
integration_rule
auto integration_rule
Definition: free_surface.cpp:187
TSPrePostProc::tsPreProc
static MoFEMErrorCode tsPreProc(TS ts)
Pre process time step.
Definition: free_surface.cpp:2872
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
ElementsAndOps
Definition: child_and_parent.cpp:18
Range
FreeSurfaceOps
Definition: FreeSurfaceOps.hpp:1
DomainEleOp
MoFEM::CoreTmp< 0 >::Initialize
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
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
save_range
auto save_range
Definition: thermo_elastic.cpp:157
get_dofs_ents_all
auto get_dofs_ents_all
get entities of dofs in the problem - used for debugging
Definition: free_surface.cpp:370
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::NaturalBC::Assembly
Assembly methods.
Definition: Natural.hpp:65
TSPrePostProc::tsCtxPtr
boost::shared_ptr< TsCtx > tsCtxPtr
Definition: free_surface.cpp:466
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
CYLINDRICAL
@ CYLINDRICAL
Definition: definitions.h:130
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MoFEM::TSMethod::ts_t
PetscReal ts_t
time
Definition: LoopMethods.hpp:166
MoFEM::AssemblyType
AssemblyType
[Storage and set boundary conditions]
Definition: FormsIntegrators.hpp:104
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
FreeSurface::setParentDofs
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.
Definition: free_surface.cpp:2806
FreeSurfaceOps::OpLhsU_dU
Lhs for U dU.
Definition: FreeSurfaceOps.hpp:473
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
mu_ave
double mu_ave
Definition: free_surface.cpp:182
FreeSurfaceOps::OpNormalConstrainRhs
Definition: FreeSurfaceOps.hpp:50
MoFEM::DMoFEMGetInterfacePtr
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
Definition: DMMoFEM.cpp:414
TSPrePostProc::tsSetIJacobian
static MoFEMErrorCode tsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Wrapper for SNES Lhs.
Definition: free_surface.cpp:3025
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
OpDomainMassU
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, U_FIELD_DIM > OpDomainMassU
Definition: free_surface.cpp:110
refine_overlap
int refine_overlap
Definition: free_surface.cpp:144
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
BITREFLEVEL_SIZE
#define BITREFLEVEL_SIZE
max number of refinements
Definition: definitions.h:232
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
FTensor::Ddg
Definition: Ddg_value.hpp:7
MoFEM::createVecScatter
auto createVecScatter(Vec x, IS ix, Vec y, IS iy)
Create a Vec Scatter object.
Definition: PetscSmartObj.hpp:361
coord_type
constexpr CoordinateTypes coord_type
Definition: incompressible_elasticity.cpp:19
MoFEM::DMMoFEMTSSetMonitor
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:1056
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
FreeSurface::solveSystem
MoFEMErrorCode solveSystem()
[Solve]
Definition: free_surface.cpp:2044
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
MoFEM::getProblemPtr
auto getProblemPtr(DM dm)
get problem pointer from DM
Definition: DMMoFEM.hpp:1044
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:64
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
get_start_bit
auto get_start_bit
Definition: free_surface.cpp:148
FreeSurface::readMesh
MoFEMErrorCode readMesh()
[Run programme]
Definition: free_surface.cpp:536
LAST_COORDINATE_SYSTEM
@ LAST_COORDINATE_SYSTEM
Definition: definitions.h:132
d_cut_off
auto d_cut_off
Definition: free_surface.cpp:206
get_M2
auto get_M2
Definition: free_surface.cpp:229
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
Monitor::postProc
boost::shared_ptr< PostProcEle > postProc
Definition: dynamic_first_order_con_law.cpp:874
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition: DeprecatedPetsc.hpp:203
FreeSurfaceOps::OpRhsU
Rhs for U.
Definition: FreeSurfaceOps.hpp:356
FreeSurfaceOps::OpCalculateLift
Definition: FreeSurfaceOps.hpp:3
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
MoFEM::MoFEMSNESMonitorFields
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
Definition: SnesCtx.cpp:239
MoFEM::FieldEntity::getHiBitNumberUId
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
Definition: FieldEntsMultiIndices.hpp:228
TSPrePostProc::globSol
SmartPetscObj< Vec > globSol
Definition: free_surface.cpp:422
TSPrePostProc::getSubVector
SmartPetscObj< Vec > getSubVector()
Definition: free_surface.cpp:3127
a0
constexpr double a0
Definition: hcurl_check_approx_in_2d.cpp:37
FreeSurfaceOps::OpLhsU_dH
Lhs for U dH.
Definition: FreeSurfaceOps.hpp:611
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
sdf_wavy_2d.w
int w
Definition: sdf_wavy_2d.py:6
Example::mField
MoFEM::Interface & mField
Definition: plastic.cpp:226
TSPrePostProc::pcSetup
static MoFEMErrorCode pcSetup(PC pc)
Definition: free_surface.cpp:3077
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
Monitor::Monitor
Monitor(SmartPetscObj< DM > dm, MoFEM::Interface &m_field, boost::shared_ptr< PostProcEle > post_proc, boost::shared_ptr< PostProcFaceEle > post_proc_bdry, boost::shared_ptr< MatrixDouble > velocity_field_ptr, boost::shared_ptr< MatrixDouble > x2_field_ptr, boost::shared_ptr< MatrixDouble > geometry_field_ptr, std::array< double, 3 > pass_field_eval_coords, boost::shared_ptr< SetPtsData > pass_field_eval_data)
Definition: dynamic_first_order_con_law.cpp:776
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MoFEM::getDMSnesCtx
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition: DMMoFEM.hpp:1127
MoFEM::FEMethod::numeredEntFiniteElementPtr
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
Definition: LoopMethods.hpp:387
marker
auto marker
set bit to marker
Definition: hanging_node_approx.cpp:82
BoundaryParentEle
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
Definition: child_and_parent.cpp:41
FreeSurface::mField
MoFEM::Interface & mField
Definition: free_surface.cpp:479
QUIET
@ QUIET
Definition: definitions.h:221
MoFEM::Unique_mi_tag
Definition: TagMultiIndices.hpp:18
mu_p
double mu_p
Definition: free_surface.cpp:166
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:238
FreeSurfaceOps::OpLhsG_dH
Lhs for H dH.
Definition: FreeSurfaceOps.hpp:1249
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
MoFEM::BasicMethod::getCacheWeakPtr
boost::weak_ptr< CacheTuple > getCacheWeakPtr() const
Get the cache weak ptr object.
Definition: LoopMethods.hpp:348
Monitor::postProcEdge
boost::shared_ptr< PostProcEleBdyCont > postProcEdge
Definition: free_surface.cpp:2038
MoFEM::SmartPetscObj< VecScatter >
FreeSurface::refineMesh
MoFEMErrorCode refineMesh(size_t overlap)
Definition: free_surface.cpp:2575
Monitor::liftVec
boost::shared_ptr< VectorDouble > liftVec
Definition: free_surface.cpp:2040
get_M0
auto get_M0
Definition: free_surface.cpp:226
FreeSurfaceOps::OpNormalConstrainLhs
Definition: FreeSurfaceOps.hpp:195
FreeSurface::setupProblem
MoFEMErrorCode setupProblem()
[Read mesh]
Definition: free_surface.cpp:552
OpDomainAssembleScalar
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpDomainAssembleScalar
Definition: free_surface.cpp:121
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEM::DMoFEMLoopFiniteElements
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:586
convert.int
int
Definition: convert.py:64
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
TSPrePostProc::globRes
SmartPetscObj< Vec > globRes
Definition: free_surface.cpp:459
MoFEM::OpAddParentEntData
Operator to project base functions from parent entity to child.
Definition: MeshProjectionDataOperators.hpp:66
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
bubble_device
auto bubble_device
Definition: free_surface.cpp:286
get_dofs_ents_by_field_name
auto get_dofs_ents_by_field_name
get entities of dofs in the problem - used for debugging
Definition: free_surface.cpp:342
MoFEM::getDMTsCtx
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1141
my_max
auto my_max
Definition: free_surface.cpp:203
DomianParentEle
PipelineManager::ElementsAndOpsByDim< FE_DIM >::DomianParentEle DomianParentEle
Definition: level_set.cpp:39
MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition: ForcesAndSourcesCore.hpp:1291
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
TSPrePostProc::fsRawPtr
Example * fsRawPtr
Definition: dynamic_first_order_con_law.cpp:397
MoFEM::OpPostProcMapInMoab::DataMapVec
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
Definition: PostProcBrokenMeshInMoabBase.hpp:711
OpDomainMassG
OpDomainMassH OpDomainMassG
Definition: free_surface.cpp:114
TSPrePostProc::pcApply
static MoFEMErrorCode pcApply(PC pc, Vec pc_f, Vec pc_x)
Definition: free_surface.cpp:3087
FreeSurface::boundaryCondition
MoFEMErrorCode boundaryCondition()
[Set up problem]
Definition: free_surface.cpp:663
tol
double tol
Definition: mesh_smoothing.cpp:26
F
@ F
Definition: free_surface.cpp:396
FR
FR
Definition: free_surface.cpp:396
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:709
Generated by Doxygen 1.8.17 and hosted at University of Glasgow