\[ \begin{aligned} \frac{\partial u(\mathbf{x}, t)}{\partial t}-\Delta u(\mathbf{x}, t) &=f(\mathbf{x}, t) & & \forall \mathbf{x} \in \Omega, t \in(0, T), \\ u(\mathbf{x}, 0) &=u_{0}(\mathbf{x}) & & \forall \mathbf{x} \in \Omega, \\ u(\mathbf{x}, t) &=g(\mathbf{x}, t) & & \forall \mathbf{x} \in \partial \Omega, t \in(0, T). \end{aligned} \]
#include <stdlib.h>
#include <cmath>
static char help[] =
"...\n\n";
PETSC>::LinearForm<GAUSS>::OpBaseTimesScalar<1>;
PETSC>::LinearForm<GAUSS>::OpSource<1, 1>;
PETSC>::LinearForm<GAUSS>::OpBaseTimesScalar<1>;
PETSC>::LinearForm<GAUSS>::OpSource<1, 1>;
: dM(dm), postProc(post_proc){};
static constexpr int saveEveryNthStep = 1;
if (ts_step % saveEveryNthStep == 0) {
"out_level_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
}
}
private:
boost::shared_ptr<PostProcEle> postProc;
};
public:
private:
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
};
}
}
};
}
if (!inner_surface.empty()) {
Range inner_surface_verts;
inner_surface_verts, false);
init_u, MBVERTEX, inner_surface_verts,
"U");
}
}
}
CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"ESSENTIAL",
"U", 0, 0);
auto &bc_map = bc_mng->getBcMapByBlockName();
for (auto b : bc_map) {
if (std::regex_match(b.first, std::regex("(.*)ESSENTIAL(.*)"))) {
for (
int i = 0;
i != b.second->bcMarkers.size(); ++
i) {
(*boundaryMarker)[
i] |= b.second->bcMarkers[
i];
}
}
}
}
auto add_domain_lhs_ops = [&](auto &pipeline) {
"U",
"U", [](
double,
double,
double) ->
double {
return k; }));
return c * fe_domain_lhs->ts_a;
};
};
auto add_domain_rhs_ops = [&](auto &pipeline) {
auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
auto dot_u_at_gauss_pts = boost::make_shared<VectorDouble>();
"U", grad_u_at_gauss_pts));
pipeline.push_back(
"U", grad_u_at_gauss_pts,
[](
double,
double,
double) ->
double {
return k; }));
"U", dot_u_at_gauss_pts,
[](
const double,
const double,
const double) {
return c; }));
auto source_term = [&](const double x, const double y, const double z) {
const auto t = fe_domain_lhs->ts_t;
return 1e1 * pow(M_E, -M_PI * M_PI *
t) * sin(1. * M_PI * x) *
sin(2. * M_PI * y);
};
};
auto add_boundary_lhs_ops = [&](auto &pipeline) {
"U",
"U", [](
const double,
const double,
const double) {
return c; }));
};
auto add_boundary_rhs_ops = [&](auto &pipeline) {
auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
auto boundary_function = [&](const double x, const double y,
const double z) {
const auto t = fe_rhs->ts_t;
return 0;
};
"U", u_at_gauss_pts,
[](
const double,
const double,
const double) {
return c; }));
};
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());
}
static PetscErrorCode
set(TS ts, PetscReal
t,
Vec u,
Vec u_t, PetscReal
a,
Mat
A, Mat B,
void *ctx) {
}
}
private:
};
auto create_post_process_element = [&]() {
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(
post_proc_fe->getOpPtrVector().push_back(
auto u_ptr = boost::make_shared<VectorDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
{{"U", u_ptr}},
{}, {}, {}
)
);
return post_proc_fe;
};
auto set_time_monitor = [&](auto dm, auto solver) {
boost::shared_ptr<Monitor> monitor_ptr(
new Monitor(dm, create_post_process_element()));
boost::shared_ptr<ForcesAndSourcesCore> null;
monitor_ptr, null, null);
};
auto set_fieldsplit_preconditioner = [&](auto solver) {
SNES snes;
CHKERR TSGetSNES(solver, &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 bc_mng = mField.getInterface<
BcManager>();
auto name_prb =
simple->getProblemName();
auto is_all_bc = bc_mng->getBlockIS(name_prb, "ESSENTIAL", "U", 0, 0);
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 set_user_ts_jacobian = [&](auto dm) {
};
auto solver = pipeline_mng->createTSIM(
CHKERR set_user_ts_jacobian(dm);
CHKERR set_time_monitor(dm, solver);
CHKERR TSSetFromOptions(solver);
CHKERR set_fieldsplit_preconditioner(solver);
}
}
}
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 heat_problem.runProgram();
}
return 0;
}