v0.14.0
shallow_wave.cpp

Solving shallow wave equation on manifold

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

/**
* \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;
#include <boost/math/quadrature/gauss_kronrod.hpp>
using namespace boost::math::quadrature;
template <int DIM> struct ElementsAndOps {};
template <> struct ElementsAndOps<2> {
};
constexpr int FE_DIM = 2;
// Use forms iterators for Grad-Grad term
// 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,
boost::shared_ptr<MatrixDouble> grad_u_ptr,
boost::shared_ptr<MatrixDouble> grad_h_ptr)
uPtr(u_ptr), uDotPtr(u_dot_ptr), uGradPtr(grad_u_ptr),
hGradPtr(grad_h_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_grad_u = getFTensor2FromMat<3, 3>(*uGradPtr);
auto t_grad_h = getFTensor1FromMat<3>(*hGradPtr);
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_grad(m, j) = t_Q(m, i) * (mu * t_grad_u(i, j));
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_grad_u;
++t_grad_h;
++t_coords;
++t_normal;
}
}
private:
boost::shared_ptr<MatrixDouble> uPtr;
boost::shared_ptr<MatrixDouble> uDotPtr;
boost::shared_ptr<MatrixDouble> uGradPtr;
boost::shared_ptr<MatrixDouble> hGradPtr;
};
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,
boost::shared_ptr<MatrixDouble> grad_u_ptr)
: AssemblyDomainEleOp(field_name_row, field_name_col,
AssemblyDomainEleOp::OPROWCOL),
uPtr(u_ptr), uGradPtr(grad_u_ptr) {
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 t_grad_u = getFTensor2FromMat<3, 3>(*uGradPtr);
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;
++t_grad_u;
}
}
private:
boost::shared_ptr<MatrixDouble> uPtr;
boost::shared_ptr<MatrixDouble> uGradPtr;
};
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 readMesh();
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 readMesh();
CHKERR setupProblem();
CHKERR createCommonData();
CHKERR boundaryCondition();
CHKERR assembleSystem();
CHKERR solveSystem();
CHKERR outputResults();
CHKERR checkResults();
}
//! [Run problem]
//! [Read mesh]
auto simple = mField.getInterface<Simple>();
CHKERR simple->getOptions();
CHKERR simple->loadFile();
}
//! [Read mesh]
//! [Set up problem]
Simple *simple = mField.getInterface<Simple>();
// Add field
CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, 3);
CHKERR simple->addDomainField("H", H1, AINSWORTH_LEGENDRE_BASE, 1);
CHKERR simple->addDataField("HO_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE, 3);
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;
PetscViewerBinaryOpen(PETSC_COMM_WORLD, "vector.dat", FILE_MODE_READ,
&viewer);
auto T = createDMVector(simple->getDM());
VecLoad(T, viewer);
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 PCFieldSplitSetType(pc, PC_COMPOSITE_ADDITIVE);
}
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 grad_u_ptr = boost::make_shared<MatrixDouble>();
auto div_u_ptr = boost::make_shared<VectorDouble>();
auto dot_h_ptr = boost::make_shared<VectorDouble>();
auto grad_h_ptr = boost::make_shared<MatrixDouble>();
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(
new OpCalculateVectorFieldGradient<3, 3>("U", grad_u_ptr));
pipeline.push_back(
pipeline.push_back(new OpCalculateScalarFieldGradient<3>("H", grad_h_ptr));
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(new OpConvectiveH("H", u_ptr, grad_h_ptr));
pipeline.push_back(
new OpURhs("U", u_ptr, dot_u_ptr, grad_u_ptr, grad_h_ptr));
};
auto set_domain_lhs = [&](auto &pipeline) {
auto u_ptr = boost::make_shared<MatrixDouble>();
auto grad_u_ptr = boost::make_shared<MatrixDouble>();
auto grad_h_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(new OpCalculateVectorFieldValues<3>("U", u_ptr));
pipeline.push_back(
new OpCalculateVectorFieldGradient<3, 3>("U", grad_u_ptr));
pipeline.push_back(new OpCalculateScalarFieldGradient<3>("H", grad_h_ptr));
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_dU("U", "U", u_ptr, grad_u_ptr));
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>();
auto grad_u_ptr = boost::make_shared<MatrixDouble>();
auto grad_h_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 OpCalculateVectorFieldGradient<3, 3>("U", grad_u_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldGradient<3>("H", grad_h_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}, {"GRAD_H", grad_h_ptr}},
{{"GRAD_U", grad_u_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();
core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "SW"));
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]
}
}
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
Example::checkResults
MoFEMErrorCode checkResults()
[Postprocess results]
Definition: dynamic_first_order_con_law.cpp:1205
g
constexpr double g
Definition: shallow_wave.cpp:63
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
FTensor::Tensor1< double, 3 >
MoFEM::OpCalculateHOJacForFaceEmbeddedIn3DSpace
OpCalculateHOJacForFaceImpl< 3 > OpCalculateHOJacForFaceEmbeddedIn3DSpace
Definition: HODataOperators.hpp:265
OpURhs
Definition: shallow_wave.cpp:82
OpBaseDivU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixScalarTimesDiv< 3 > OpBaseDivU
Definition: shallow_wave.cpp:36
MoFEM::OpCalculateDivergenceVectorFieldValues
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:486
EntData
EntitiesFieldData::EntData EntData
Definition: shallow_wave.cpp:26
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
init_h
auto init_h
Initialisation function.
Definition: free_surface.cpp:295
help
static char help[]
[Check]
Definition: shallow_wave.cpp:900
phi
static double phi
Definition: poisson_2d_dis_galerkin.cpp:30
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
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::OpCalculateScalarFieldValuesDot
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
Definition: UserDataOperators.hpp:273
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
i
FTensor::Index< 'i', 3 > i
Definition: shallow_wave.cpp:77
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM.hpp
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
OpConvectiveH_dGradH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::TriLinearForm< GAUSS >::OpConvectiveTermLhsDy< 1, 1, 3 > OpConvectiveH_dGradH
Definition: shallow_wave.cpp:60
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
OpBase
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
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::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
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1293
ROW
@ ROW
Definition: definitions.h:136
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
OpBaseGradH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixVectorTimesGrad< 1, 3, 3 > OpBaseGradH
Definition: shallow_wave.cpp:38
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:1502
Example::readMesh
MoFEMErrorCode readMesh()
[Run problem]
Definition: dynamic_first_order_con_law.cpp:463
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:667
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:548
u_max
constexpr double u_max
Definition: shallow_wave.cpp:68
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
main
int main(int argc, char *argv[])
Definition: shallow_wave.cpp:902
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
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
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:317
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::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
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
j
FTensor::Index< 'j', 3 > j
Definition: shallow_wave.cpp:78
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
MoFEM::OpGetHONormalsOnFace
Calculate normals at Gauss points of triangle element.
Definition: HODataOperators.hpp:281
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
beta_montain
constexpr double beta_montain
Definition: shallow_wave.cpp:73
penalty
constexpr double penalty
Definition: shallow_wave.cpp:75
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
OpBaseTimesDivU
OpBaseTimesDotH OpBaseTimesDivU
Definition: shallow_wave.cpp:48
BiLinearForm
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
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1535
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:385
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
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
Monitor
[Push operators to pipeline]
Definition: dynamic_first_order_con_law.cpp:783
approx_order
int approx_order
Definition: test_broken_space.cpp:50
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
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::OpInvertMatrix
Definition: UserDataOperators.hpp:3249
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
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
h_hat
constexpr double h_hat
Definition: shallow_wave.cpp:67
omega
constexpr double omega
Definition: shallow_wave.cpp:62
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
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::SmartPetscObj< IS >
l
FTensor::Index< 'l', 3 > l
Definition: shallow_wave.cpp:79
FE_DIM
constexpr int FE_DIM
Definition: shallow_wave.cpp:24
Example::outputResults
MoFEMErrorCode outputResults()
[Solve]
Definition: dynamic_first_order_con_law.cpp:1183
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
MoFEM::OpSetInvJacH1ForFaceEmbeddedIn3DSpace
Definition: UserDataOperators.hpp:2932
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
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:359
OpBaseDivU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixScalarTimesDiv< SPACE_DIM > OpBaseDivU
[Essential boundary conditions]
Definition: seepage.cpp:79
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
mu
constexpr double mu
Definition: shallow_wave.cpp:64
OpBaseTimesDotH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBaseTimesDotH
Definition: shallow_wave.cpp:47