#ifndef __NONLINEARPOISSON2D_HPP__
#define __NONLINEARPOISSON2D_HPP__
#include <stdlib.h>
PETSC>::LinearForm<GAUSS>::OpBaseTimesScalar<1>;
PETSC>::LinearForm<GAUSS>::OpSource<1, 1>;
typedef boost::function<
double(
const double,
const double,
const double)>
public:
OpDomainLhs(
std::string row_field_name, std::string col_field_name,
boost::shared_ptr<VectorDouble> field_vec,
boost::shared_ptr<MatrixDouble> field_grad_mat)
fieldVec(field_vec), fieldGradMat(field_grad_mat) {}
const int nb_row_dofs = row_data.
getIndices().size();
const int nb_col_dofs = col_data.
getIndices().size();
const double area = getMeasure();
const int nb_integration_points = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
auto t_field_grad = getFTensor1FromMat<2>(*fieldGradMat);
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * area;
for (int rr = 0; rr != nb_row_dofs; ++rr) {
for (int cc = 0; cc != nb_col_dofs; cc++) {
locLhs(rr, cc) += (((1 + t_field * t_field) * t_row_diff_base(
i) *
(2.0 * t_field * t_field_grad(
i) *
t_row_diff_base(
i) * t_col_base)) *
++t_col_base;
++t_col_diff_base;
}
++t_row_base;
++t_row_diff_base;
}
++t_w;
++t_field;
++t_field_grad;
}
}
private:
boost::shared_ptr<VectorDouble> fieldVec;
boost::shared_ptr<MatrixDouble> fieldGradMat;
};
public:
OpDomainRhs(
boost::shared_ptr<VectorDouble> field_vec,
boost::shared_ptr<MatrixDouble> field_grad_mat)
sourceTermFunc(source_term_function), fieldVec(field_vec),
fieldGradMat(field_grad_mat) {}
const double area = getMeasure();
const int nb_integration_points = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
auto t_coords = getFTensor1CoordsAtGaussPts();
auto t_field_grad = getFTensor1FromMat<2>(*fieldGradMat);
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * area;
double body_source =
sourceTermFunc(t_coords(0), t_coords(1), t_coords(2));
nf[rr] +=
(-t_base * body_source +
t_grad_base(
i) * t_field_grad(
i) * (1 + t_field * t_field)) *
++t_base;
++t_grad_base;
}
++t_w;
++t_coords;
++t_field;
++t_field_grad;
}
}
private:
boost::shared_ptr<VectorDouble> fieldVec;
boost::shared_ptr<MatrixDouble> fieldGradMat;
};
};
#endif //__NONLINEARPOISSON2D_HPP__
#include <stdlib.h>
static char help[] =
"...\n\n";
inline double sqr(
const double x) {
return x * x; }
inline double cube(
const double x) {
return x * x * x; }
public:
private:
static double sourceTermFunction(const double x, const double y,
const double z) {
return 2 * M_PI * M_PI *
(cos(M_PI * x) * cos(M_PI * y) +
cube(cos(M_PI * x)) *
cube(cos(M_PI * y)) -
cos(M_PI * x) * cos(M_PI * y) *
(
sqr(sin(M_PI * x)) *
sqr(cos(M_PI * y)) +
sqr(sin(M_PI * y)) *
sqr(cos(M_PI * x))));
}
static double boundaryFunction(const double x, const double y,
const double z) {
return -cos(M_PI * x) *
cos(M_PI * y);
}
int atomTest = 0;
enum NORMS { NORM = 0, LAST_NORM };
enum EX_NORMS { EX_NORM = 0, LAST_EX_NORM };
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
boost::shared_ptr<PostProcEle> postProc;
Range boundaryEntitiesForFieldsplit;
};
: mField(m_field) {}
}
}
true);
auto get_ents_by_dim = [&](const auto dim) {
return domain_ents;
} else {
if (dim == 0)
else
return ents;
}
};
auto get_base = [&]() {
auto domain_ents = get_ents_by_dim(
SPACE_DIM);
if (domain_ents.empty())
case MBQUAD:
case MBHEX:
case MBTRI:
case MBTET:
default:
}
};
auto base = get_base();
PETSC_NULL);
}
auto domain_rule_lhs = [](
int,
int,
int p) ->
int {
return 2 * p - 1; };
auto domain_rule_rhs = [](
int,
int,
int p) ->
int {
return 2 * p - 1; };
auto boundary_rule_lhs = [](
int,
int,
int p) ->
int {
return 2 * p; };
auto boundary_rule_rhs = [](
int,
int,
int p) ->
int {
return 2 * p; };
CHKERR pipeline_mng->setDomainRhsIntegrationRule(domain_rule_lhs);
CHKERR pipeline_mng->setDomainLhsIntegrationRule(domain_rule_rhs);
CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(boundary_rule_lhs);
CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(boundary_rule_rhs);
}
auto get_ents_on_mesh_skin = [&]() {
std::string entity_name = it->getName();
if (entity_name.compare(0, 18, "BOUNDARY_CONDITION") == 0) {
boundary_entities, true);
}
}
boundary_vertices, true);
boundary_entities.merge(boundary_vertices);
return boundary_entities;
};
auto mark_boundary_dofs = [&](
Range &&skin_edges) {
auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
ProblemsManager::OR, skin_edges, *marker_ptr);
return marker_ptr;
};
}
auto add_domain_lhs_ops = [&](auto &pipeline) {
auto data_u_at_gauss_pts = boost::make_shared<VectorDouble>();
auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
pipeline.push_back(
pipeline.push_back(new OpDomainLhs(
};
auto add_domain_rhs_ops = [&](auto &pipeline) {
auto data_u_at_gauss_pts = boost::make_shared<VectorDouble>();
auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
pipeline.push_back(
data_u_at_gauss_pts,
grad_u_at_gauss_pts));
};
auto add_boundary_lhs_ops = [&](auto &pipeline) {
[](const double, const double, const double) { return 1; }));
};
auto add_boundary_rhs_ops = [&](auto &pipeline) {
auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
pipeline.push_back(
[](const double, const double, const double) { return 1; }));
};
add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
add_boundary_lhs_ops(pipeline_mng->getOpBoundaryLhsPipeline());
add_boundary_rhs_ops(pipeline_mng->getOpBoundaryRhsPipeline());
}
auto set_fieldsplit_preconditioner = [&](auto snes) {
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
PC pc;
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs == PETSC_TRUE) {
auto name_prb =
simple->getProblemName();
int is_all_bc_size;
CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
<< "Field split block size " << is_all_bc_size;
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
is_all_bc);
}
};
auto solver = pipeline_mng->createSNES();
CHKERR SNESSetFromOptions(solver);
CHKERR set_fieldsplit_preconditioner(solver);
CHKERR SNESSolve(solver, global_rhs, global_solution);
SCATTER_REVERSE);
}
auto post_proc = boost::make_shared<PostProcEle>(
mField);
auto u_ptr = boost::make_shared<VectorDouble>();
post_proc->getOpPtrVector().push_back(
post_proc->getOpPtrVector().push_back(
new OpPPMap(post_proc->getPostProcMesh(), post_proc->getMapGaussPts(),
{{domainField, u_ptr}},
{}, {}, {}
)
);
CHKERR post_proc->writeFile(
"out_result.h5m");
}
auto check_result_fe_ptr = boost::make_shared<DomainEle>(
mField);
check_result_fe_ptr->getOpPtrVector(), {H1})),
"Apply transform");
check_result_fe_ptr->getRuleHook = [](
int,
int,
int p) {
return p; };
auto analyticalFunction = [&](double x, double y, double z) {
return cos(M_PI * x) * cos(M_PI * y);
};
auto u_ptr = boost::make_shared<VectorDouble>();
check_result_fe_ptr->getOpPtrVector().push_back(
auto mValFuncPtr = boost::make_shared<VectorDouble>();
check_result_fe_ptr->getOpPtrVector().push_back(
check_result_fe_ptr->getOpPtrVector().push_back(
check_result_fe_ptr->getOpPtrVector().push_back(
check_result_fe_ptr);
const double *L2_norms;
const double *Ex_norms;
<<
"L2_NORM: " << std::sqrt(L2_norms[
NORM]);
<< "NORMALISED_ERROR: "
<< (std::sqrt(L2_norms[
NORM]) / std::sqrt(Ex_norms[
EX_NORM]));
}
const double *L2_norms;
const double *Ex_norms;
double ref_percentage = 4.4e-04;
double normalised_error;
case 1:
normalised_error = std::sqrt(L2_norms[0]) / std::sqrt(Ex_norms[0]);
break;
default:
"atom test %d does not exist",
atomTest);
}
if (normalised_error > ref_percentage) {
"atom test %d failed! Calculated Norm %3.16e is greater than "
"reference Norm %3.16e",
atomTest, normalised_error, ref_percentage);
}
}
}
int main(
int argc,
char *argv[]) {
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");
try {
DMType dm_name = "DMMOFEM";
CHKERR poisson_problem.runProgram();
}
return 0;
}