v0.13.0
VEC-4: Shallow wave equation on sphere
Note
Prerequisites of this tutorial include VEC-0: Linear elasticity


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

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 [43].

/**
* \file shallow_wave.cpp
* \example shallow_wave.cpp
*
* Solving shallow wave equation on manifold
*
* The inital conditions are set following this paper \cite scott2016test.
*
*/
/* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
#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;
FormsIntegrators<DomainEleOp>::Assembly<PETSC>::OpBase;
// Use forms iterators for Grad-Grad term
using OpBaseDivU = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
using OpBaseGradH = FormsIntegrators<DomainEleOp>::Assembly<
// Use forms for Mass term
using OpMassUU = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
using OpMassHH = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
using OpBaseTimesDotH = FormsIntegrators<DomainEleOp>::Assembly<
using OpSourceU = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
GAUSS>::OpSource<1, 3>;
using OpSourceH = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
GAUSS>::OpSource<1, 1>;
using OpConvectiveH = FormsIntegrators<DomainEleOp>::Assembly<
using OpConvectiveH_dU = FormsIntegrators<DomainEleOp>::Assembly<
using OpConvectiveH_dGradH = FormsIntegrators<DomainEleOp>::Assembly<
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)
: AssemblyDomainEleOp(field_name, field_name, AssemblyDomainEleOp::OPROW),
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 = smartCreateDMVector(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);
post_proc_fe->generateReferenceElementMesh();
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->addFieldValuesPostProc("U");
post_proc_fe->addFieldValuesPostProc("H");
post_proc_fe->addFieldValuesPostProc("HO_POSITIONS");
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();
SmartPetscObj<IS> is_u;
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 = smartCreateDMVector(dm);
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(
new OpCalculateDivergenceVectorFieldValues<3>("U", div_u_ptr));
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:
SmartPetscObj<DM> dM;
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();
SmartPetscObj<IS> is_u;
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);
post_proc_fe->generateReferenceElementMesh();
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(
new OpSetInvJacH1ForFaceEmbeddedIn3DSpace(inv_jac_ptr));
post_proc_fe->addFieldValuesPostProc("U");
post_proc_fe->addFieldValuesPostProc("H");
post_proc_fe->addFieldValuesGradientPostProc("U");
post_proc_fe->addFieldValuesGradientPostProc("H");
post_proc_fe->addFieldValuesPostProc("HO_POSITIONS");
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 = smartCreateDMVector(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"));
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
int main(int argc, char *argv[])
static char help[]
constexpr double a
EntitiesFieldData::EntData EntData
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Kronecker Delta class.
Tensor1< T, Tensor_Dim > normalize()
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
constexpr double omega
constexpr double alpha
@ ROW
Definition: definitions.h:136
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ H1
continuous field
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
const double init_u
auto integration_rule
constexpr double lambda
auto init_h
constexpr auto t_kd
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixScalarTimesDiv< SPACE_DIM, coord_type > OpMixScalarTimesDiv
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:956
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:481
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:544
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:364
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:311
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:342
FTensor::Index< 'i', SPACE_DIM > i
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:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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: DMMMoFEM.cpp:994
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)
CoreTmp< 0 > Core
Definition: Core.hpp:1096
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
OpCalculateHOJacForFaceImpl< 3 > OpCalculateHOJacForFaceEmbeddedIn3DSpace
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
const double D
diffusivity
const double r
rate factor
const double u0
inital vale on blocksets
int save_every_nth_step
static double phi
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:41
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 penalty
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpBaseTimesDotH
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
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpConvectiveTermLhsDy< 1, 1, 3 > OpConvectiveH_dGradH
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]
MoFEMErrorCode boundaryCondition()
[Set up problem]
MoFEMErrorCode assembleSystem()
[Boundary condition]
MoFEMErrorCode readMesh()
[Run problem]
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode solveSystem()
[Solve]
MoFEMErrorCode createCommonData()
[Create common data]
MoFEMErrorCode runProblem()
[Create common data]
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode outputResults()
[Solve]
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:279
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:323
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
[Push operators to pipeline]
Postprocess on face.

Higher order geometry

/**
* \file approx_sphere.cpp
* \example approx_sphere.cpp
*
*/
/* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
#include <BasicFiniteElements.hpp>
constexpr int FM_DIM = 2;
template <int DIM> struct ElementsAndOps {};
template <> struct ElementsAndOps<2> {
};
FormsIntegrators<DomainEleOp>::Assembly<PETSC>::OpBase;
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)
: AssemblyDomainEleOp(field_name, field_name, AssemblyDomainEleOp::OPROW),
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(field_name, field_name,
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)
: DomainEleOp(field_name, field_name, AssemblyDomainEleOp::OPROW),
xPtr(x_ptr) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
auto t_w = getFTensor0IntegrationWeight();
auto t_x = getFTensor1FromMat<3>(*xPtr);
auto t_normal = getFTensor1NormalsAtGaussPts();
auto nb_integration_pts = getGaussPts().size2();
double error = 0;
for (int gg = 0; gg != nb_integration_pts; gg++) {
double l = std::sqrt(t_normal(i) * t_normal(i));
error += t_w * l * std::abs((t_x(i) * t_x(i) - A * A));
++t_w;
++t_x;
++t_normal;
}
CHKERR VecSetValue(errorVec, 0, error, ADD_VALUES);
}
static SmartPetscObj<Vec> errorVec;
private:
boost::shared_ptr<MatrixDouble> xPtr;
};
SmartPetscObj<Vec> OpError::errorVec;
struct ApproxSphere {
ApproxSphere(MoFEM::Interface &m_field) : mField(m_field) {}
MoFEMErrorCode runProblem();
private:
MoFEMErrorCode getOptions();
MoFEMErrorCode readMesh();
MoFEMErrorCode setupProblem();
MoFEMErrorCode setOPs();
MoFEMErrorCode solveSystem();
MoFEMErrorCode outputResults();
};
//! [Run programme]
CHKERR getOptions();
CHKERR readMesh();
CHKERR setupProblem();
CHKERR setOPs();
CHKERR solveSystem();
CHKERR outputResults();
}
//! [Run programme]
}
//! [Read mesh]
auto simple = mField.getInterface<Simple>();
CHKERR simple->getOptions();
CHKERR simple->loadFile();
}
//! [Read mesh]
//! [Set up problem]
auto simple = mField.getInterface<Simple>();
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]
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
auto integration_rule = [](int, int, int approx_order) {
return 3 * approx_order;
};
CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
auto x_ptr = boost::make_shared<MatrixDouble>();
auto dot_x_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto def_ops = [&](auto &pipeline) {
pipeline.push_back(
new OpCalculateVectorFieldValues<3>("HO_POSITIONS", x_ptr));
pipeline.push_back(
new OpCalculateVectorFieldValuesDot<3>("HO_POSITIONS", dot_x_ptr));
};
def_ops(pipeline_mng->getOpDomainRhsPipeline());
def_ops(pipeline_mng->getOpDomainLhsPipeline());
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpRhs("HO_POSITIONS", x_ptr, dot_x_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpLhs("HO_POSITIONS", x_ptr, dot_x_ptr));
}
//! [Push operators to pipeline]
//! [Solve]
// Project HO geometry from mesh
Projection10NodeCoordsOnField ent_method_material(mField, "HO_POSITIONS");
CHKERR mField.loop_dofs("HO_POSITIONS", ent_method_material);
auto *simple = mField.getInterface<Simple>();
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto dm = simple->getDM();
ts = pipeline_mng->createTSIM();
double ftime = 1;
CHKERR TSSetMaxSteps(ts, 1);
CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
auto T = smartCreateDMVector(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 simple = mField.getInterface<Simple>();
auto dm = simple->getDM();
auto post_proc_fe = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
post_proc_fe->generateReferenceElementMesh();
post_proc_fe->addFieldValuesPostProc("HO_POSITIONS");
CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
CHKERR post_proc_fe->writeFile("out_approx.h5m");
auto error_fe = boost::make_shared<DomainEle>(mField);
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>();
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::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]
}
}
auto lhs_J_dx2
constexpr int FM_DIM
auto res_J_dx
auto res_J
constexpr double a2
constexpr double a4
FTensor::Index< 'k', 3 > k
auto createSmartVectorMPI
Create MPI Vector.
double A
MoFEMErrorCode readMesh()
[Read mesh]
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode runProblem()
[Run programme]
MoFEMErrorCode setOPs()
[Set up problem]
MoFEMErrorCode getOptions()
[Run programme]
MoFEMErrorCode outputResults()
[Solve]
static SmartPetscObj< Vec > errorVec