v0.14.0
Loading...
Searching...
No Matches
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

Example solution

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 [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
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixScalarTimesDiv< SPACE_DIM > OpBaseDivU
Definition: seepage.cpp:80
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 >::TriLinearForm< GAUSS >::OpConvectiveTermLhsDy< 1, 1, 3 > OpConvectiveH_dGradH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 3 > OpMassUU
constexpr double h0
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpConvectiveTermRhs< 1, 1, 3 > OpConvectiveH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::TriLinearForm< GAUSS >::OpConvectiveTermLhsDu< 1, 1, 3 > OpConvectiveH_dU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixScalarTimesDiv< 3 > OpBaseDivU
OpBaseTimesDotH OpBaseTimesDivU
constexpr double u_max
FTensor::Index< 'j', 3 > j
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:226
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:440
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:262
MoFEM::Interface & mField
Definition: plastic.cpp:233
boost::shared_ptr< FEMethod > domianRhsFEPtr
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:274
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

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_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;
}
CHKERR VecSetValue(errorVec, 0, error, ADD_VALUES);
}
private:
boost::shared_ptr<MatrixDouble> xPtr;
};
struct ApproxSphere {
ApproxSphere(MoFEM::Interface &m_field) : mField(m_field) {}
private:
};
//! [Run programme]
}
//! [Run programme]
}
//! [Read mesh]
CHKERR simple->getOptions();
CHKERR simple->loadFile();
}
//! [Read mesh]
//! [Set up problem]
CHKERR simple->addDomainField("HO_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE, 3);
CHKERR simple->addDataField("HO_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE, 3);
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]
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 *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 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";
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
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]
}
}
constexpr double A
auto lhs_J_dx2
constexpr int FM_DIM
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
auto res_J_dx
auto res_J
FTensor::Index< 'k', 3 > k
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
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.
constexpr double a2
constexpr double a4
FTensor::Index< 'k', 3 > k
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
constexpr AssemblyType A
MoFEMErrorCode readMesh()
[Read mesh]
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
MoFEM::Interface & mField
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode runProblem()
[Run programme]
MoFEMErrorCode setOPs()
[Set up problem]
MoFEMErrorCode getOptions()
[Run programme]
MoFEMErrorCode outputResults()
[Solve]
virtual int get_comm_size() const =0
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Basic algebra on fields.
Definition: FieldBlas.hpp:21
auto getFTensor0IntegrationWeight()
Get integration weights.
@ OPROW
operator doWork function is executed on FE rows
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Projection of edge entities with one mid-node on hierarchical basis.
static SmartPetscObj< Vec > errorVec
boost::shared_ptr< MatrixDouble > xPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate grad-grad operator.
boost::shared_ptr< MatrixDouble > xPtr
boost::shared_ptr< MatrixDouble > xDotPtr
boost::shared_ptr< MatrixDouble > xPtr
boost::shared_ptr< MatrixDouble > xDotPtr
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data)
Class dedicated to integrate operator.