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

Solving shallow wave equation on manifold

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

/**
* \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
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) {}
private:
boost::shared_ptr<FEMethod> domianLhsFEPtr;
boost::shared_ptr<FEMethod> domianRhsFEPtr;
};
//! [Create common data]
}
//! [Create common data]
//! [Run problem]
}
//! [Run problem]
//! [Read mesh]
CHKERR simple->getOptions();
CHKERR simple->loadFile();
}
//! [Read mesh]
//! [Set up problem]
// 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 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 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 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){};
constexpr int save_every_nth_step = 50;
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 *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";
// 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]
}
}
std::string param_file
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
static char help[]
int main()
Definition: adol-c_atom.cpp:46
constexpr double a
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Kronecker Delta class.
Tensor1< T, Tensor_Dim > normalize()
@ ROW
Definition: definitions.h:123
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
constexpr int order
constexpr double omega
Save field DOFS on vertices/tags.
const double init_u
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
static double penalty
FTensor::Index< 'm', SPACE_DIM > m
static double phi
@ F
auto integration_rule
auto init_h
Initialisation function.
constexpr auto t_kd
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMixScalarTimesDiv< SPACE_DIM, COORD_TYPE > OpMixScalarTimesDiv
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:509
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:572
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1003
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
FTensor::Index< 'i', SPACE_DIM > i
static double lambda
double D
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
const double T
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition: DMMoFEM.cpp:1042
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
const double u0
inital vale on blocksets
int r
Definition: sdf.py:8
int save_every_nth_step
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr auto field_name
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
constexpr int FE_DIM
constexpr double mu
constexpr double phi_0
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpMassHH
constexpr double phi_1
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 3 > OpSourceU
constexpr double beta_montain
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 3 > OpMassUU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpConvectiveTermLhsDu< 1, 1, 3 > OpConvectiveH_dU
constexpr double h0
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpConvectiveTermRhs< 1, 1, 3 > OpConvectiveH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixScalarTimesDiv< 3 > OpBaseDivU
OpBaseTimesDotH OpBaseTimesDivU
constexpr double u_max
FTensor::Index< 'j', 3 > j
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpConvectiveTermLhsDy< 1, 1, 3 > OpConvectiveH_dGradH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBaseTimesDotH
constexpr double omega
FTensor::Index< 'i', 3 > i
constexpr double g
FTensor::Index< 'm', 3 > m
constexpr double h_hat
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpSourceH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixVectorTimesGrad< 1, 3, 3 > OpBaseGradH
constexpr double phi_2
constexpr double alpha_montain
[Example]
Definition: plastic.cpp:222
MoFEMErrorCode boundaryCondition()
[Set up problem]
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
MoFEMErrorCode readMesh()
[Run problem]
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode solveSystem()
[Solve]
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:428
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:258
MoFEM::Interface & mField
Definition: plastic.cpp:229
boost::shared_ptr< FEMethod > domianRhsFEPtr
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:270
MoFEMErrorCode outputResults()
[Solve]
boost::shared_ptr< FEMethod > domianLhsFEPtr
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
VectorDouble locF
local entity vector
int nbRows
number of dofs on rows
int nbIntegrationPts
number of integration points
MatrixDouble locMat
local entity block matrix
int nbCols
number if dof on column
int nbRowBaseFunctions
number or row base functions
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
Calculate HO coordinates at gauss points.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get rate of scalar field at integration points.
Get value at integration points for scalar field.
Approximate field values for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Calculate normals at Gauss points of triangle element.
Post post-proc data at points from hash maps.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
intrusive_ptr for managing petsc objects
PetscInt ts_step
time step number
Vec & ts_u
state vector
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
[Push operators to pipeline]
SmartPetscObj< DM > dM
MoFEMErrorCode postProcess()
function is run at the end of loop
boost::shared_ptr< PostProcEle > postProc
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate grad-grad operator.
boost::shared_ptr< MatrixDouble > uPtr
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate grad-grad operator.
boost::shared_ptr< MatrixDouble > uGradPtr
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data)
Class dedicated to integrate operator.
boost::shared_ptr< MatrixDouble > uDotPtr
boost::shared_ptr< MatrixDouble > uGradPtr
boost::shared_ptr< MatrixDouble > hGradPtr
boost::shared_ptr< MatrixDouble > uPtr