v0.14.0
VEC-4: Shallow wave equation on sphere

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

Note
Intended learning outcome:
• solving PhD on mainfold
• handling two vectorial fields
• setting inital boundary conditions
• material time derivative
• hyperbolic equation
• approximation sphere with high-order polynomials

# Supplement lectures

You can find many resources on the internet about shallow wave equations. We point to two of them, which has been used as a starting point for this tutorial.

# Plain program

The inital conditions are set following this paper [63].

/**
* \file shallow_wave.cpp
* \example shallow_wave.cpp
*
* Solving shallow wave equation on manifold
*
* The inital conditions are set following this paper \cite scott2016test.
*
*/
#include <MoFEM.hpp>
using namespace MoFEM;
template <int DIM> struct ElementsAndOps {};
template <> struct ElementsAndOps<2> {
};
constexpr int FE_DIM = 2;
// Use forms for Mass term
PETSC>::LinearForm<GAUSS>::OpBaseTimesScalar<1>;
GAUSS>::OpSource<1, 3>;
GAUSS>::OpSource<1, 1>;
constexpr double omega = 7.292 * 1e-5;
constexpr double g = 9.80616;
constexpr double mu = 1e4;
constexpr double h0 = 1e4;
constexpr double h_hat = 120;
constexpr double u_max = 80;
constexpr double phi_0 = M_PI / 7;
constexpr double phi_1 = M_PI / 2 - phi_0;
constexpr double phi_2 = M_PI / 4;
constexpr double alpha_montain = 1. / 3.;
constexpr double beta_montain = 1. / 15.;
constexpr double penalty = 1;
struct OpURhs : public AssemblyDomainEleOp {
OpURhs(const std::string field_name, boost::shared_ptr<MatrixDouble> u_ptr,
boost::shared_ptr<MatrixDouble> u_dot_ptr,
const double vol = getMeasure();
auto t_w = getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor0N();
auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
auto t_dot_u = getFTensor1FromMat<3>(*uDotPtr);
auto t_u = getFTensor1FromMat<3>(*uPtr);
auto t_coords = getFTensor1CoordsAtGaussPts();
auto t_normal = getFTensor1NormalsAtGaussPts();
for (int gg = 0; gg != nbIntegrationPts; gg++) {
const double alpha = t_w * vol;
auto t_nf = getFTensor1FromArray<3, 3>(locF);
const auto a = std::sqrt(t_coords(i) * t_coords(i));
const auto sin_fi = t_coords(2) / a;
const auto f = 2 * omega * sin_fi;
t_r(i) = t_normal(i);
t_r.normalize();
t_P(i, j) = t_r(i) * t_r(j);
t_Q(i, j) = t_kd(i, j) - t_P(i, j);
t_A(i, m) = levi_civita(i, j, m) * t_r(j);
t_rhs(m) = t_Q(m, i) * (t_dot_u(i) + t_grad_u(i, j) * t_u(j) +
f * t_A(i, j) * t_u(j) + g * t_grad_h(i));
t_rhs(m) += t_P(m, j) * t_u(j);
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
t_nf(i) += alpha * t_row_base * t_rhs(i);
t_nf(i) += alpha * t_row_diff_base(j) * t_rhs_grad(i, j);
++t_row_base;
++t_row_diff_base;
++t_nf;
}
for (; rr < nbRowBaseFunctions; ++rr) {
++t_row_base;
++t_row_diff_base;
}
++t_w;
++t_u;
++t_dot_u;
++t_coords;
++t_normal;
}
}
private:
boost::shared_ptr<MatrixDouble> uPtr;
boost::shared_ptr<MatrixDouble> uDotPtr;
};
struct OpULhs_dU : public AssemblyDomainEleOp {
OpULhs_dU(const std::string field_name_row, const std::string field_name_col,
boost::shared_ptr<MatrixDouble> u_ptr,
: AssemblyDomainEleOp(field_name_row, field_name_col,
AssemblyDomainEleOp::OPROWCOL),
this->sYmm = false;
}
const double vol = getMeasure();
auto t_w = getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor0N();
auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
auto t_coords = getFTensor1CoordsAtGaussPts();
auto t_normal = getFTensor1NormalsAtGaussPts();
auto t_u = getFTensor1FromMat<3>(*uPtr);
auto get_t_mat = [&](const int rr) {
&locMat(rr + 0, 0), &locMat(rr + 0, 1), &locMat(rr + 0, 2),
&locMat(rr + 1, 0), &locMat(rr + 1, 1), &locMat(rr + 1, 2),
&locMat(rr + 2, 0), &locMat(rr + 2, 1), &locMat(rr + 2, 2)};
};
const auto ts_a = getFEMethod()->ts_a;
for (int gg = 0; gg != nbIntegrationPts; gg++) {
const auto a = std::sqrt(t_coords(i) * t_coords(i));
const auto sin_fi = t_coords(2) / a;
const auto f = 2 * omega * sin_fi;
t_r(i) = t_normal(i);
t_r.normalize();
t_P(i, j) = t_r(i) * t_r(j);
t_Q(i, j) = t_kd(i, j) - t_P(i, j);
t_A(i, m) = levi_civita(i, j, m) * t_r(j);
t_rhs_du(m, j) =
t_Q(m, i) * (ts_a * t_kd(i, j) + t_grad_u(i, j) + f * t_A(i, j)) +
t_P(m, j);
const double alpha = t_w * vol;
int rr = 0;
for (; rr != nbRows / 3; rr++) {
auto t_col_base = col_data.getFTensor0N(gg, 0);
auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
auto t_mat = get_t_mat(3 * rr);
for (int cc = 0; cc != nbCols / 3; cc++) {
t_mat(i, j) += (alpha * t_row_base * t_col_base) * t_rhs_du(i, j);
t_mat(i, j) += (alpha * mu) * t_Q(i, j) *
(t_row_diff_base(m) * t_col_diff_base(m));
++t_col_diff_base;
++t_col_base;
++t_mat;
}
++t_row_base;
++t_row_diff_base;
}
for (; rr < nbRowBaseFunctions; ++rr) {
++t_row_base;
++t_row_diff_base;
}
++t_w;
++t_coords;
++t_normal;
++t_u;
}
}
private:
boost::shared_ptr<MatrixDouble> uPtr;
};
struct OpULhs_dH : public AssemblyDomainEleOp {
OpULhs_dH(const std::string field_name_row, const std::string field_name_col)
: AssemblyDomainEleOp(field_name_row, field_name_col,
AssemblyDomainEleOp::OPROWCOL) {
this->sYmm = false;
}
// get element volume
const double vol = getMeasure();
// get integration weights
auto t_w = getFTensor0IntegrationWeight();
// get base function gradient on rows
auto t_row_base = row_data.getFTensor0N();
// normal
auto t_normal = getFTensor1NormalsAtGaussPts();
auto get_t_vec = [&](const int rr) {
&locMat(rr + 0, 0),
&locMat(rr + 1, 0),
&locMat(rr + 2, 0)};
};
// loop over integration points
for (int gg = 0; gg != nbIntegrationPts; gg++) {
t_r(i) = t_normal(i);
t_r.normalize();
t_P(i, j) = t_r(i) * t_r(j);
t_Q(i, j) = t_kd(i, j) - t_P(i, j);
const double alpha = t_w * vol;
int rr = 0;
for (; rr != nbRows / 3; rr++) {
auto t_vec = get_t_vec(3 * rr);
auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
const double a = alpha * g * t_row_base;
for (int cc = 0; cc != nbCols; cc++) {
t_vec(i) += a * (t_Q(i, m) * t_col_diff_base(m));
++t_vec;
++t_col_diff_base;
}
++t_row_base;
}
for (; rr < nbRowBaseFunctions; ++rr)
++t_row_base;
++t_w;
++t_normal;
}
}
};
struct Example {
Example(MoFEM::Interface &m_field) : mField(m_field) {}
MoFEMErrorCode runProblem();
private:
MoFEMErrorCode setupProblem();
MoFEMErrorCode createCommonData();
MoFEMErrorCode boundaryCondition();
MoFEMErrorCode assembleSystem();
MoFEMErrorCode solveSystem();
MoFEMErrorCode outputResults();
MoFEMErrorCode checkResults();
boost::shared_ptr<FEMethod> domianLhsFEPtr;
boost::shared_ptr<FEMethod> domianRhsFEPtr;
};
//! [Create common data]
}
//! [Create common data]
//! [Run problem]
CHKERR setupProblem();
CHKERR createCommonData();
CHKERR boundaryCondition();
CHKERR assembleSystem();
CHKERR solveSystem();
CHKERR outputResults();
CHKERR checkResults();
}
//! [Run problem]
auto simple = mField.getInterface<Simple>();
CHKERR simple->getOptions();
}
//! [Set up problem]
int order = 2;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
CHKERR simple->setFieldOrder("U", order);
CHKERR simple->setFieldOrder("H", order);
CHKERR simple->setFieldOrder("HO_POSITIONS", 3);
CHKERR simple->setUp();
}
//! [Set up problem]
//! [Boundary condition]
PetscBool is_restart = PETSC_FALSE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_restart", &is_restart,
PETSC_NULL);
auto restart_vector = [&]() {
auto simple = mField.getInterface<Simple>();
auto dm = simple->getDM();
MOFEM_LOG("SW", Sev::inform)
<< "reading vector in binary from vector.dat ...";
PetscViewer viewer;
&viewer);
auto T = createDMVector(simple->getDM());
CHKERR DMoFEMMeshToLocalVector(dm, T, INSERT_VALUES, SCATTER_REVERSE);
};
if (is_restart) {
CHKERR restart_vector();
} else {
const double e_n = exp(-4 / pow(phi_1 - phi_0, 2));
const double u0 = u_max / e_n;
FTensor::Tensor1<double, 3> t_k{0., 0., 1.};
t_A(i, m) = levi_civita(i, j, m) * t_k(j);
auto get_phi = [&](const double x, const double y, const double z) {
const double r = std::sqrt(t_r(i) * t_r(i));
return atan2(z, r);
};
auto init_u_phi = [&](const double phi) {
if (phi > phi_0 && phi < phi_1) {
return u0 * exp(1. / ((phi - phi_0) * (phi - phi_1)));
} else {
return 0.;
}
};
auto init_u = [&](const double x, const double y, const double z) {
FTensor::Tensor1<double, 3> t_u{0., 0., 0.};
const double u_phi = init_u_phi(get_phi(x, y, z));
if (u_phi > 0) {
t_r.normalize();
t_u(i) = ((t_A(i, j) * t_r(j)) * u_phi);
}
return t_u;
};
auto init_h = [&](const double x, const double y, const double z) {
const double a = std::sqrt(x * x + y * y + z * z);
auto integral = [&](const double fi) {
const double u_phi = init_u_phi(fi);
const auto f = 2 * omega * sin(fi);
return a * u_phi * (f + (tan(fi) / a) * u_phi);
};
auto montain = [&](const double lambda, const double fi) {
if (lambda > -M_PI && lambda < M_PI)
return h_hat * cos(fi) * exp(-pow(lambda / alpha_montain, 2)) *
exp(-pow((phi_2 - fi) / beta_montain, 2));
else
return 0.;
};
const double fi = get_phi(x, y, z);
const double lambda = atan2(x, y);
double h1 = 0;
if (fi > phi_0)
h1 = gauss_kronrod<double, 32>::integrate(
integral, phi_0, fi, 0, std::numeric_limits<float>::epsilon());
return h0 + (montain(lambda, fi) - (h1 / g));
};
auto set_domain_general = [&](auto &pipeline) {
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(new OpGetHONormalsOnFace("HO_POSITIONS"));
pipeline.push_back(new OpCalculateHOCoords("HO_POSITIONS"));
pipeline.push_back(new OpSetHOWeightsOnFace());
};
auto set_domain_rhs = [&](auto &pipeline) {
pipeline.push_back(new OpSourceU("U", init_u));
pipeline.push_back(new OpSourceH("H", init_h));
};
auto set_domain_lhs = [&](auto &pipeline) {
pipeline.push_back(
new OpMassUU("U", "U", [](double, double, double) { return 1; }));
pipeline.push_back(
new OpMassHH("H", "H", [](double, double, double) { return 1; }));
};
auto post_proc = [&]() {
auto simple = mField.getInterface<Simple>();
auto dm = simple->getDM();
auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
new OpGetHONormalsOnFace("HO_POSITIONS"));
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateHOCoords("HO_POSITIONS"));
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
auto u_ptr = boost::make_shared<MatrixDouble>();
auto h_ptr = boost::make_shared<VectorDouble>();
auto pos_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("H", h_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>("HO_POSITIONS", pos_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
{{"H", h_ptr}},
{{"U", u_ptr}, {"HO_POSITIONS", pos_ptr}},
{}, {}
)
);
CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
CHKERR post_proc_fe->writeFile("out_init.h5m");
};
auto solve_init = [&]() {
auto simple = mField.getInterface<Simple>();
auto pipeline_mng = mField.getInterface<PipelineManager>();
auto solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(solver);
PC pc;
CHKERR KSPGetPC(solver, &pc);
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs == PETSC_TRUE) {
auto bc_mng = mField.getInterface<BcManager>();
auto name_prb = simple->getProblemName();
CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
name_prb, ROW, "U", 0, 3, is_u);
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_u);
}
CHKERR KSPSetUp(solver);
auto dm = simple->getDM();
auto D = createDMVector(dm);
auto F = vectorDuplicate(D);
CHKERR KSPSolve(solver, F, D);
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
};
auto pipeline_mng = mField.getInterface<PipelineManager>();
auto integration_rule = [](int, int, int approx_order) {
return 2 * approx_order + 4;
};
CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
set_domain_general(pipeline_mng->getOpDomainRhsPipeline());
set_domain_rhs(pipeline_mng->getOpDomainRhsPipeline());
set_domain_general(pipeline_mng->getOpDomainLhsPipeline());
set_domain_lhs(pipeline_mng->getOpDomainLhsPipeline());
CHKERR solve_init();
CHKERR post_proc();
// Clear pipelines
pipeline_mng->getOpDomainRhsPipeline().clear();
pipeline_mng->getOpDomainLhsPipeline().clear();
}
}
//! [Boundary condition]
//! [Push operators to pipeline]
// Push element from reference configuration to current configuration in 3d
// space
auto set_domain_general = [&](auto &pipeline) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(new OpGetHONormalsOnFace("HO_POSITIONS"));
pipeline.push_back(new OpCalculateHOCoords("HO_POSITIONS"));
pipeline.push_back(new OpSetHOWeightsOnFace());
pipeline.push_back(new OpCalculateHOJacForFaceEmbeddedIn3DSpace(jac_ptr));
pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(new OpSetInvJacH1ForFaceEmbeddedIn3DSpace(inv_jac_ptr));
};
auto set_domain_rhs = [&](auto &pipeline) {
auto dot_u_ptr = boost::make_shared<MatrixDouble>();
auto u_ptr = boost::make_shared<MatrixDouble>();
auto div_u_ptr = boost::make_shared<VectorDouble>();
auto dot_h_ptr = boost::make_shared<VectorDouble>();
pipeline.push_back(new OpCalculateVectorFieldValuesDot<3>("U", dot_u_ptr));
pipeline.push_back(new OpCalculateScalarFieldValuesDot("H", dot_h_ptr));
pipeline.push_back(new OpCalculateVectorFieldValues<3>("U", u_ptr));
pipeline.push_back(
pipeline.push_back(
pipeline.push_back(new OpBaseTimesDotH(
"H", dot_h_ptr, [](double, double, double) { return 1.; }));
pipeline.push_back(new OpBaseTimesDivU(
"H", div_u_ptr, [](double, double, double) { return h0; }));
pipeline.push_back(
};
auto set_domain_lhs = [&](auto &pipeline) {
auto u_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(new OpCalculateVectorFieldValues<3>("U", u_ptr));
pipeline.push_back(
pipeline.push_back(new OpMassHH("H", "H", [&](double, double, double) {
return domianLhsFEPtr->ts_a;
}));
pipeline.push_back(new OpBaseDivU(
"H", "U", [](const double, const double, const double) { return h0; },
false, false));
pipeline.push_back(
new OpConvectiveH_dU("H", "U", grad_h_ptr, []() { return 1; }));
pipeline.push_back(
new OpConvectiveH_dGradH("H", "H", u_ptr, []() { return 1; }));
pipeline.push_back(new OpULhs_dH("U", "H"));
};
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto integration_rule = [](int, int, int approx_order) {
return 2 * approx_order + 4;
};
CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
set_domain_general(pipeline_mng->getOpDomainRhsPipeline());
set_domain_general(pipeline_mng->getOpDomainLhsPipeline());
set_domain_rhs(pipeline_mng->getOpDomainRhsPipeline());
set_domain_lhs(pipeline_mng->getOpDomainLhsPipeline());
domianLhsFEPtr = pipeline_mng->getDomainLhsFE();
domianRhsFEPtr = pipeline_mng->getDomainRhsFE();
}
//! [Push operators to pipeline]
/**
* @brief Monitor solution
*
* This functions is called by TS solver at the end of each step. It is used
* to output results to the hard drive.
*/
struct Monitor : public FEMethod {
Monitor(SmartPetscObj<DM> dm, boost::shared_ptr<PostProcEle> post_proc)
: dM(dm), postProc(post_proc){};
MoFEMErrorCode postProcess() {
constexpr int save_every_nth_step = 50;
if (ts_step % save_every_nth_step == 0) {
CHKERR DMoFEMLoopFiniteElements(dM, "dFE", postProc);
CHKERR postProc->writeFile(
"out_step_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
MOFEM_LOG("SW", Sev::verbose)
<< "writing vector in binary to vector.dat ...";
PetscViewer viewer;
PetscViewerBinaryOpen(PETSC_COMM_WORLD, "vector.dat", FILE_MODE_WRITE,
&viewer);
VecView(ts_u, viewer);
PetscViewerDestroy(&viewer);
}
}
private:
boost::shared_ptr<PostProcEle> postProc;
};
//! [Solve]
auto *simple = mField.getInterface<Simple>();
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto dm = simple->getDM();
auto set_initial_step = [&](auto ts) {
int step = 0;
CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-step", &step,
PETSC_NULL);
CHKERR TSSetStepNumber(ts, step);
};
auto set_fieldsplit_preconditioner_ksp = [&](auto ksp) {
PC pc;
CHKERR KSPGetPC(ksp, &pc);
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs == PETSC_TRUE) {
auto bc_mng = mField.getInterface<BcManager>();
auto name_prb = simple->getProblemName();
CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
name_prb, ROW, "U", 0, 3, is_u);
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_u);
}
};
// Setup postprocessing
auto get_fe_post_proc = [&]() {
auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
new OpGetHONormalsOnFace("HO_POSITIONS"));
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateHOCoords("HO_POSITIONS"));
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
post_proc_fe->getOpPtrVector().push_back(
auto u_ptr = boost::make_shared<MatrixDouble>();
auto h_ptr = boost::make_shared<VectorDouble>();
auto pos_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("H", h_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>("HO_POSITIONS", pos_ptr));
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"H", h_ptr}},
)
);
return post_proc_fe;
};
auto set_fieldsplit_preconditioner_ts = [&](auto solver) {
SNES snes;
CHKERR TSGetSNES(solver, &snes);
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
CHKERR set_fieldsplit_preconditioner_ksp(ksp);
};
ts = pipeline_mng->createTSIM();
boost::shared_ptr<FEMethod> null_fe;
auto monitor_ptr = boost::make_shared<Monitor>(dm, get_fe_post_proc());
CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
null_fe, monitor_ptr);
// Add monitor to time solver
double ftime = 1;
// CHKERR TSSetMaxTime(ts, ftime);
CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
auto T = createDMVector(simple->getDM());
CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR TSSetSolution(ts, T);
CHKERR TSSetFromOptions(ts);
CHKERR set_fieldsplit_preconditioner_ts(ts);
CHKERR TSSetUp(ts);
CHKERR set_initial_step(ts);
CHKERR TSSolve(ts, NULL);
CHKERR TSGetTime(ts, &ftime);
}
//! [Solve]
//! [Postprocess results]
}
//! [Postprocess results]
//! [Check]
}
//! [Check]
static char help[] = "...\n\n";
int main(int argc, char *argv[]) {
// 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();
LogManager::setLog("SW");
MOFEM_LOG_TAG("SW", "example");
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
//! [Register MoFEM discrete manager in PETSc
//! [Create MoAB]
moab::Core mb_instance; ///< mesh database
moab::Interface &moab = mb_instance; ///< mesh database interface
//! [Create MoAB]
//! [Create MoFEM]
MoFEM::Core core(moab); ///< finite element database
MoFEM::Interface &m_field = core; ///< finite element database insterface
//! [Create MoFEM]
//! [Example]
Example ex(m_field);
CHKERR ex.runProblem();
//! [Example]
}
}

# Higher order geometry

/**
* \file approx_sphere.cpp
* \example approx_sphere.cpp
*
*/
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
constexpr int FM_DIM = 2;
template <int DIM> struct ElementsAndOps {};
template <> struct ElementsAndOps<2> {
};
constexpr double a = 1;
constexpr double a2 = a * a;
constexpr double a4 = a2 * a2;
constexpr double A = 6371220;
auto res_J = [](const double x, const double y, const double z) {
const double res = (x * x + y * y + z * z - a2);
return res;
};
auto res_J_dx = [](const double x, const double y, const double z) {
const double res = res_J(x, y, z);
return FTensor::Tensor1<double, 3>{res * (2 * x), res * (2 * y),
res * (2 * z)};
};
auto lhs_J_dx2 = [](const double x, const double y, const double z) {
const double res = res_J(x, y, z);
(res * 2 + (4 * x * x)),
(4 * y * x),
(4 * z * x),
(4 * x * y),
(2 * res + (4 * y * y)),
(4 * z * y),
(4 * x * z),
(4 * y * z),
(2 * res + (4 * z * z))};
};
struct OpRhs : public AssemblyDomainEleOp {
OpRhs(const std::string field_name, boost::shared_ptr<MatrixDouble> x_ptr,
boost::shared_ptr<MatrixDouble> dot_x_ptr)
xPtr(x_ptr), xDotPtr(dot_x_ptr) {}
auto t_w = getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor0N();
auto t_x0 = getFTensor1CoordsAtGaussPts();
auto t_x = getFTensor1FromMat<3>(*xPtr);
auto t_dot_x = getFTensor1FromMat<3>(*xDotPtr);
auto t_normal = getFTensor1NormalsAtGaussPts();
for (int gg = 0; gg != nbIntegrationPts; gg++) {
FTensor::Tensor1<double, 3> t_n{t_x0(0), t_x0(1), t_x0(2)};
t_n.normalize();
t_P(i, j) = t_n(i) * t_n(j);
t_Q(i, j) = t_kd(i, j) - t_P(i, j);
auto t_J_res = res_J_dx(t_x(0), t_x(1), t_x(2));
const double alpha = t_w;
auto t_nf = getFTensor1FromArray<3, 3>(locF);
double l = std::sqrt(t_normal(i) * t_normal(i));
t_res(i) =
alpha * l * ((t_P(i, k) * t_J_res(k) + t_Q(i, k) * t_dot_x(k)));
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
t_nf(j) += t_row_base * t_res(j);
++t_row_base;
++t_nf;
}
for (; rr < nbRowBaseFunctions; ++rr) {
++t_row_base;
}
++t_w;
++t_x;
++t_dot_x;
++t_x0;
++t_normal;
}
}
private:
boost::shared_ptr<MatrixDouble> xPtr;
boost::shared_ptr<MatrixDouble> xDotPtr;
};
struct OpLhs : public AssemblyDomainEleOp {
OpLhs(const std::string field_name, boost::shared_ptr<MatrixDouble> x_ptr,
boost::shared_ptr<MatrixDouble> dot_x_ptr)
AssemblyDomainEleOp::OPROWCOL),
xPtr(x_ptr), xDotPtr(dot_x_ptr) {
this->sYmm = false;
}
auto t_w = getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor0N();
auto t_x0 = getFTensor1CoordsAtGaussPts();
auto t_x = getFTensor1FromMat<3>(*xPtr);
auto t_dot_x = getFTensor1FromMat<3>(*xDotPtr);
auto t_normal = getFTensor1NormalsAtGaussPts();
auto get_t_mat = [&](const int rr) {
&locMat(rr + 0, 0), &locMat(rr + 0, 1), &locMat(rr + 0, 2),
&locMat(rr + 1, 0), &locMat(rr + 1, 1), &locMat(rr + 1, 2),
&locMat(rr + 2, 0), &locMat(rr + 2, 1), &locMat(rr + 2, 2)};
};
const double ts_a = getTSa();
for (int gg = 0; gg != nbIntegrationPts; gg++) {
FTensor::Tensor1<double, 3> t_n{t_x0(0), t_x0(1), t_x0(2)};
t_n.normalize();
t_P(i, j) = t_n(i) * t_n(j);
t_Q(i, j) = t_kd(i, j) - t_P(i, j);
auto t_J_lhs = lhs_J_dx2(t_x(0), t_x(1), t_x(2));
double l = std::sqrt(t_normal(i) * t_normal(i));
const double alpha = t_w;
t_lhs(i, j) =
(alpha * l) * (t_P(i, k) * t_J_lhs(k, j) + t_Q(i, j) * ts_a);
int rr = 0;
for (; rr != nbRows / 3; rr++) {
auto t_col_base = col_data.getFTensor0N(gg, 0);
auto t_mat = get_t_mat(3 * rr);
for (int cc = 0; cc != nbCols / 3; cc++) {
const double rc = t_row_base * t_col_base;
t_mat(i, j) += rc * t_lhs(i, j);
++t_col_base;
++t_mat;
}
++t_row_base;
}
for (; rr < nbRowBaseFunctions; ++rr)
++t_row_base;
++t_w;
++t_x;
++t_x0;
++t_normal;
}
}
private:
boost::shared_ptr<MatrixDouble> xPtr;
boost::shared_ptr<MatrixDouble> xDotPtr;
};
struct OpError : public DomainEleOp {
OpError(const std::string field_name, boost::shared_ptr<MatrixDouble> x_ptr)
xPtr(x_ptr) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
auto t_w = getFTensor0IntegrationWeight();
auto t_x = getFTensor1FromMat<3>(*xPtr);
auto t_normal = getFTensor1NormalsAtGaussPts();
auto nb_integration_pts = getGaussPts().size2();
double error = 0;
for (int gg = 0; gg != nb_integration_pts; gg++) {
double l = std::sqrt(t_normal(i) * t_normal(i));
error += t_w * l * std::abs((t_x(i) * t_x(i) - A * A));
++t_w;
++t_x;
++t_normal;
}
}
static SmartPetscObj<Vec> errorVec;
private:
boost::shared_ptr<MatrixDouble> xPtr;
};
struct ApproxSphere {
ApproxSphere(MoFEM::Interface &m_field) : mField(m_field) {}
MoFEMErrorCode runProblem();
private:
MoFEMErrorCode getOptions();
MoFEMErrorCode setupProblem();
MoFEMErrorCode setOPs();
MoFEMErrorCode solveSystem();
MoFEMErrorCode outputResults();
};
//! [Run programme]
CHKERR getOptions();
CHKERR setupProblem();
CHKERR setOPs();
CHKERR solveSystem();
CHKERR outputResults();
}
//! [Run programme]
}
auto simple = mField.getInterface<Simple>();
CHKERR simple->getOptions();
}
//! [Set up problem]
auto simple = mField.getInterface<Simple>();
int order = 3;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
CHKERR simple->setFieldOrder("HO_POSITIONS", order);
CHKERR simple->setUp();
}
//! [Set up problem]
//! [Push operators to pipeline]
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
auto integration_rule = [](int, int, int approx_order) {
return 3 * approx_order;
};
auto x_ptr = boost::make_shared<MatrixDouble>();
auto dot_x_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto def_ops = [&](auto &pipeline) {
pipeline.push_back(
new OpCalculateVectorFieldValues<3>("HO_POSITIONS", x_ptr));
pipeline.push_back(
new OpCalculateVectorFieldValuesDot<3>("HO_POSITIONS", dot_x_ptr));
};
def_ops(pipeline_mng->getOpDomainRhsPipeline());
def_ops(pipeline_mng->getOpDomainLhsPipeline());
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpRhs("HO_POSITIONS", x_ptr, dot_x_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpLhs("HO_POSITIONS", x_ptr, dot_x_ptr));
}
//! [Push operators to pipeline]
//! [Solve]
// Project HO geometry from mesh
Projection10NodeCoordsOnField ent_method_material(mField, "HO_POSITIONS");
CHKERR mField.loop_dofs("HO_POSITIONS", ent_method_material);
auto *simple = mField.getInterface<Simple>();
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto dm = simple->getDM();
ts = pipeline_mng->createTSIM();
double ftime = 1;
CHKERR TSSetMaxSteps(ts, 1);
CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
auto T = createDMVector(simple->getDM());
CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR TSSetSolution(ts, T);
CHKERR TSSetFromOptions(ts);
CHKERR TSSolve(ts, NULL);
CHKERR TSGetTime(ts, &ftime);
CHKERR mField.getInterface<FieldBlas>()->fieldScale(A, "HO_POSITIONS");
}
//! [Solve]
auto x_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto simple = mField.getInterface<Simple>();
auto dm = simple->getDM();
auto post_proc_fe =
boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>("HO_POSITIONS", x_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
{},
{{"HO_POSITIONS", x_ptr}},
{}, {}
)
);
CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
CHKERR post_proc_fe->writeFile("out_approx.h5m");
auto error_fe = boost::make_shared<DomainEle>(mField);
error_fe->getOpPtrVector().push_back(
new OpGetHONormalsOnFace("HO_POSITIONS"));
error_fe->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>("HO_POSITIONS", x_ptr));
error_fe->getOpPtrVector().push_back(new OpError("HO_POSITIONS", x_ptr));
error_fe->preProcessHook = [&]() {
MOFEM_LOG("EXAMPLE", Sev::inform) << "Create vec ";
mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
VecZeroEntries(OpError::errorVec);
};
error_fe->postProcessHook = [&]() {
CHKERR VecAssemblyBegin(OpError::errorVec);
CHKERR VecAssemblyEnd(OpError::errorVec);
double error2;
CHKERR VecSum(OpError::errorVec, &error2);
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "Error " << std::sqrt(error2 / (4 * M_PI * A * A));
};
CHKERR DMoFEMLoopFiniteElements(dm, "dFE", error_fe);
CHKERR simple->deleteDM();
CHKERR simple->deleteFiniteElements();
if (mField.get_comm_size() > 1)
CHKERR mField.get_moab().write_file("out_ho_mesh.h5m", "MOAB",
"PARALLEL=WRITE_PART");
else
CHKERR mField.get_moab().write_file("out_ho_mesh.h5m");
}
//! [Postprocess results]
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
MoFEM::Core::Initialize(&argc, &argv, param_file, help);
auto core_log = logging::core::get();
LogManager::setLog("EXAMPLE");
MOFEM_LOG_TAG("EXAMPLE", "example");
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
//! [Register MoFEM discrete manager in PETSc
//! [Create MoAB]
moab::Core mb_instance; ///< mesh database
moab::Interface &moab = mb_instance; ///< mesh database interface
//! [Create MoAB]
//! [Create MoFEM]
MoFEM::Core core(moab); ///< finite element database
MoFEM::Interface &m_field = core; ///< finite element database insterface
//! [Create MoFEM]
//! [ApproxSphere]
ApproxSphere ex(m_field);
CHKERR ex.runProblem();
//! [ApproxSphere]
}
}
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
a2
constexpr double a2
Definition: hcurl_check_approx_in_2d.cpp:39
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
Example::checkResults
MoFEMErrorCode checkResults()
[Postprocess results]
Definition: dynamic_first_order_con_law.cpp:1205
g
constexpr double g
Definition: shallow_wave.cpp:63
omega
constexpr double omega
Save field DOFS on vertices/tags.
Definition: dynamic_first_order_con_law.cpp:93
MoFEM::CoreInterface::loop_dofs
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
Example::assembleSystem
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
Definition: dynamic_first_order_con_law.cpp:647
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
OpError
Definition: initial_diffusion.cpp:61
FTensor::Tensor1< double, 3 >
ApproxSphere::getOptions
MoFEMErrorCode getOptions()
[Run programme]
Definition: approx_sphere.cpp:292
MoFEM::OpCalculateHOJacForFaceEmbeddedIn3DSpace
OpCalculateHOJacForFaceImpl< 3 > OpCalculateHOJacForFaceEmbeddedIn3DSpace
Definition: HODataOperators.hpp:265
OpURhs
Definition: shallow_wave.cpp:82
MoFEM::OpCalculateDivergenceVectorFieldValues
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:486
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
init_h
auto init_h
Initialisation function.
Definition: free_surface.cpp:295
ApproxSphere::setupProblem
MoFEMErrorCode setupProblem()
Definition: approx_sphere.cpp:308
phi
static double phi
Definition: poisson_2d_dis_galerkin.cpp:29
help
static char help[]
Definition: activate_deactivate_dofs.cpp:13
Definition: approx_sphere.cpp:298
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
OpConvectiveTermRhs
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpConvectiveTermRhs< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermRhs
Definition: operators_tests.cpp:52
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
MoFEM::PipelineManager::getOpDomainRhsPipeline
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
Definition: PipelineManager.hpp:773
OpConvectiveTermLhsDu
FormsIntegrators< DomainEleOp >::Assembly< A >::TriLinearForm< I >::OpConvectiveTermLhsDu< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermLhsDu
Definition: operators_tests.cpp:56
phi_1
constexpr double phi_1
Definition: shallow_wave.cpp:70
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::OpCalculateScalarFieldValuesDot
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
Definition: UserDataOperators.hpp:273
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
OpError::errorVec
static SmartPetscObj< Vec > errorVec
Definition: approx_sphere.cpp:254
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM.hpp
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
OpSourceU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 3 > OpSourceU
Definition: shallow_wave.cpp:51
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
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::TriLinearForm< GAUSS >::OpConvectiveTermLhsDy< 1, 1, 3 > OpConvectiveH_dGradH
Definition: shallow_wave.cpp:60
BasicFiniteElements.hpp
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
ApproxSphere::runProblem
MoFEMErrorCode runProblem()
[Run programme]
Definition: approx_sphere.cpp:280
OpBase
OpBaseImpl< PETSC, EdgeEleOp > OpBase
FTensor::levi_civita
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
Definition: Levi_Civita.hpp:617
MoFEM::LogManager::createSink
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::OpCalculateHOCoords
Calculate HO coordinates at gauss points.
Definition: HODataOperators.hpp:41
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2< double, 3, 3 >
order
constexpr int order
Definition: dg_projection.cpp:18
OpMixScalarTimesDiv
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMixScalarTimesDiv< SPACE_DIM, COORD_TYPE > OpMixScalarTimesDiv
Definition: free_surface.cpp:125
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
h0
constexpr double h0
Definition: shallow_wave.cpp:65
OpULhs_dH
Definition: shallow_wave.cpp:267
FE_DIM
constexpr int FE_DIM
[Define dimension]
Definition: level_set.cpp:19
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1294
ROW
@ ROW
Definition: definitions.h:123
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixVectorTimesGrad< 1, 3, 3 > OpBaseGradH
Definition: shallow_wave.cpp:38
approx_order
static constexpr int approx_order
Definition: prism_polynomial_approximation.cpp:14
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl
Approximate field values for given petsc vector.
Definition: UserDataOperators.hpp:595
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1489
res_J_dx
auto res_J_dx
Definition: approx_sphere.cpp:45
[Run problem]
Definition: dynamic_first_order_con_law.cpp:463
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:178
save_every_nth_step
int save_every_nth_step
Definition: photon_diffusion.cpp:67
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
u_max
constexpr double u_max
Definition: shallow_wave.cpp:68
FM_DIM
constexpr int FM_DIM
Definition: approx_sphere.cpp:15
ApproxSphere::setOPs
MoFEMErrorCode setOPs()
[Set up problem]
Definition: approx_sphere.cpp:325
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1067
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
ApproxSphere::outputResults
MoFEMErrorCode outputResults()
[Solve]
Definition: approx_sphere.cpp:394
ReactionDiffusionEquation::u0
const double u0
inital vale on blocksets
Definition: reaction_diffusion.cpp:24
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
phi_2
constexpr double phi_2
Definition: shallow_wave.cpp:71
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
OpLhs
Definition: approx_sphere.cpp:134
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
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
Example
[Example]
Definition: plastic.cpp:177
convert.type
type
Definition: convert.py:64
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:310
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
OpConvectiveTermLhsDy
FormsIntegrators< DomainEleOp >::Assembly< A >::TriLinearForm< I >::OpConvectiveTermLhsDy< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermLhsDy
Definition: operators_tests.cpp:60
Example::solveSystem
MoFEMErrorCode solveSystem()
[Solve]
Definition: dynamic_first_order_con_law.cpp:893
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::CoreInterface::get_comm_size
virtual int get_comm_size() const =0
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
MoFEM::OpSetHOWeightsOnFace
Modify integration weights on face to take in account higher-order geometry.
Definition: HODataOperators.hpp:122
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
MoFEM::OpGetHONormalsOnFace
Calculate normals at Gauss points of triangle element.
Definition: HODataOperators.hpp:281
MoFEM::PipelineManager::setDomainRhsIntegrationRule
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:530
OpSourceH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpSourceH
Definition: shallow_wave.cpp:53
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
MoFEM::PipelineManager::getOpDomainLhsPipeline
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
Definition: PipelineManager.hpp:749
beta_montain
constexpr double beta_montain
Definition: shallow_wave.cpp:73
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
penalty
constexpr double penalty
Definition: shallow_wave.cpp:75
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
lhs_J_dx2
auto lhs_J_dx2
Definition: approx_sphere.cpp:51
a4
constexpr double a4
Definition: hcurl_check_approx_in_2d.cpp:41
ApproxSphere::solveSystem
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
Definition: approx_sphere.cpp:361
OpBaseTimesDivU
OpBaseTimesDotH OpBaseTimesDivU
Definition: shallow_wave.cpp:48
BiLinearForm
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', 3 >
FTensor::Tensor1::normalize
Tensor1< T, Tensor_Dim > normalize()
Definition: Tensor1_value.hpp:77
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1536
OpRhs
Definition: approx_sphere.cpp:68
MoFEM::PipelineManager::setDomainLhsIntegrationRule
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:503
integration_rule
auto integration_rule
Definition: free_surface.cpp:185
ElementsAndOps
Definition: child_and_parent.cpp:18
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_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::EntitiesFieldData::EntData::getFTensor1DiffN
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition: EntitiesFieldData.cpp:526
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1102
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
Example::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:225
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
Example::boundaryCondition
MoFEMErrorCode boundaryCondition()
[Set up problem]
Definition: dynamic_first_order_con_law.cpp:536
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
Monitor
[Push operators to pipeline]
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
mu
double mu
Definition: dynamic_first_order_con_law.cpp:98
OpConvectiveH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpConvectiveTermRhs< 1, 1, 3 > OpConvectiveH
Definition: shallow_wave.cpp:56
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
OpConvectiveH_dU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::TriLinearForm< GAUSS >::OpConvectiveTermLhsDu< 1, 1, 3 > OpConvectiveH_dU
Definition: shallow_wave.cpp:58
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition: PetscSmartObj.hpp:202
MoFEM::OpInvertMatrix
Definition: UserDataOperators.hpp:3254
Example::runProblem
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:213
init_u
constexpr double init_u
Definition: heat_equation.cpp:58
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
OpMassHH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpMassHH
Definition: shallow_wave.cpp:44
ApproxSphere
Definition: approx_sphere.cpp:262
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
h_hat
constexpr double h_hat
Definition: shallow_wave.cpp:67
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
alpha_montain
constexpr double alpha_montain
Definition: shallow_wave.cpp:72
Example::createCommonData
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:396
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
MoFEM::SmartPetscObj< IS >
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
Example::outputResults
MoFEMErrorCode outputResults()
[Solve]
Definition: dynamic_first_order_con_law.cpp:1183
res_J
auto res_J
Definition: approx_sphere.cpp:40
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:416
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
OpMassUU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 3 > OpMassUU
Definition: shallow_wave.cpp:42
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
OpBaseDivU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixScalarTimesDiv< SPACE_DIM > OpBaseDivU
Definition: seepage.cpp:80
phi_0
constexpr double phi_0
Definition: shallow_wave.cpp:69
OpULhs_dU
Definition: shallow_wave.cpp:168
F
@ F
Definition: free_surface.cpp:394
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
OpBaseTimesDotH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBaseTimesDotH
Definition: shallow_wave.cpp:47