TODO: Add more operators.
static char help[] =
"...\n\n";
};
};
template <int FIELD_DIM>
template <int FIELD_DIM>
template <int FIELD_DIM>
template <int FIELD_DIM>
template <int FIELD_DIM>
constexpr
bool debug =
false;
int main(
int argc,
char *argv[]) {
try {
DMType dm_name = "DMMOFEM";
auto core_log = logging::core::get();
core_log->add_sink(
auto post_proc = [&](auto dm, auto f_res, auto out_name) {
auto post_proc_fe =
boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(m_field);
auto scal_vec = boost::make_shared<VectorDouble>();
auto vec_mat = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
f_res));
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"SCALAR", scal_vec}},
{{"VECTOR", vec_mat}},
{}, {})
);
post_proc_fe);
post_proc_fe->writeFile(out_name);
};
auto TestOpGradGrad = [&]() {
MOFEM_LOG(
"OpTester", Sev::verbose) <<
"TestOpGradGrad";
pip->getOpDomainLhsPipeline().clear();
pip->getOpDomainRhsPipeline().clear();
pip->setDomainLhsIntegrationRule([](int, int, int o) { return 2 * o; });
pip->setDomainLhsIntegrationRule([](int, int, int o) { return 2 * o; });
"SCALAR", "SCALAR", [](double, double, double) { return 1; }));
"VECTOR", "VECTOR", [](double, double, double) { return 1; }));
auto scl_mat = boost::make_shared<MatrixDouble>();
auto vec_mat = boost::make_shared<MatrixDouble>();
pip->getOpDomainRhsPipeline().push_back(
pip->getOpDomainRhsPipeline().push_back(
vec_mat));
"SCALAR", scl_mat, [](double, double, double) { return 1; }));
"VECTOR", vec_mat, [](double, double, double) { return 1; }));
constexpr
double eps = 1e-6;
auto x = opt->setRandomFields(
simple->getDM(),
{{"SCALAR", {-1, 1}}, {"VECTOR", {-1, 1}}});
auto diff_x = opt->setRandomFields(
simple->getDM(), {{
"SCALAR", {-1, 1}}, {
"VECTOR", {-1, 1}}});
auto diff_res = opt->checkCentralFiniteDifference(
simple->getDM(),
simple->getDomainFEName(), pip->getDomainRhsFE(),
double fnorm;
CHKERR VecNorm(diff_res, NORM_2, &fnorm);
MOFEM_LOG_C(
"OpTester", Sev::inform,
"TestOpGradGrad %3.4e", fnorm);
constexpr double err = 1e-9;
if (fnorm > err)
"Norm of directional derivative too large");
opt->directionalCentralFiniteDifference(
"out_directional_directive_gradgrad.h5m");
}
};
auto TestOpConvection = [&]() {
MOFEM_LOG(
"OpTester", Sev::verbose) <<
"TestOpConvection";
pip->getOpDomainLhsPipeline().clear();
pip->getOpDomainRhsPipeline().clear();
pip->setDomainLhsIntegrationRule([](int, int, int o) { return 2 * o; });
pip->setDomainLhsIntegrationRule([](int, int, int o) { return 2 * o; });
auto evaluate_field = [&](auto &pipe) {
auto scl_mat = boost::make_shared<MatrixDouble>();
auto vec_mat = boost::make_shared<MatrixDouble>();
auto dot_vec_vel = boost::make_shared<MatrixDouble>();
pipe.push_back(
"VECTOR", vec_mat));
"VECTOR", dot_vec_vel));
return std::make_tuple(scl_mat, vec_mat, dot_vec_vel);
};
auto [rhs_scl_mat, rhs_vec_mat, rhs_dot_vec_vel] =
evaluate_field(pip->getOpDomainRhsPipeline());
pip->getOpDomainRhsPipeline().push_back(
pip->getOpDomainRhsPipeline().push_back(
rhs_vec_mat));
auto [lhs_scl_mat, lhs_vec_mat, lhs_dot_vec_vel] =
evaluate_field(pip->getOpDomainLhsPipeline());
auto op_convective_term_lhs_du_scalar =
op_convective_term_lhs_du_scalar->feScalingFun =
[](
const FEMethod *fe_ptr) {
return fe_ptr->ts_a; };
pip->getOpDomainLhsPipeline().push_back(op_convective_term_lhs_du_scalar);
pip->getOpDomainLhsPipeline().push_back(
auto op_convective_term_lhs_du_vec =
op_convective_term_lhs_du_vec->feScalingFun = [](
const FEMethod *fe_ptr) {
return fe_ptr->ts_a;
};
pip->getOpDomainLhsPipeline().push_back(op_convective_term_lhs_du_vec);
pip->getOpDomainLhsPipeline().push_back(
lhs_dot_vec_vel));
constexpr
double eps = 1e-6;
auto x = opt->setRandomFields(
simple->getDM(),
{{"SCALAR", {-1, 1}}, {"VECTOR", {-1, 1}}});
auto x_t = opt->setRandomFields(
simple->getDM(), {{
"SCALAR", {-1, 1}}, {
"VECTOR", {-1, 1}}});
auto diff_x = opt->setRandomFields(
simple->getDM(), {{
"SCALAR", {-1, 1}}, {
"VECTOR", {-1, 1}}});
auto diff_res = opt->checkCentralFiniteDifference(
simple->getDM(),
simple->getDomainFEName(), pip->getDomainRhsFE(),
double fnorm;
CHKERR VecNorm(diff_res, NORM_2, &fnorm);
MOFEM_LOG_C(
"OpTester", Sev::inform,
"TestOpConvection %3.4e", fnorm);
constexpr double err = 1e-9;
if (fnorm > err) {
"Norm of directional derivative too large");
}
opt->directionalCentralFiniteDifference(
pip->getDomainRhsFE(), x, x_t,
"out_directional_directive_convection.h5m");
}
};
}
return 0;
}