static char help[] =
"...\n\n";
static double fUn(const double x, const double y, double z) {
for (
int i = 0;
i <= o; ++
i) {
r += pow(x,
i) * pow(y,
j);
}
}
}
double z) {
for (
int i = 0;
i <= o; ++
i) {
r(0) +=
i > 0 ?
i * pow(x,
i - 1) * pow(y,
j) : 0;
r(1) +=
j > 0 ?
j * pow(x,
i) * pow(y,
j - 1) : 0;
}
}
}
}
};
static double fUn(const double x, const double y, double z) {
for (
int i = 0;
i <= o; ++
i) {
for (
int j = 0;
j <= o -
i;
j++) {
r += pow(x,
i) * pow(y,
j) * pow(z,
k);
}
}
}
}
}
double z) {
for (
int i = 0;
i <= o; ++
i) {
for (
int j = 0;
j <= o -
i;
j++) {
r(0) +=
i > 0 ?
i * pow(x,
i - 1) * pow(y,
j) * pow(z,
k) : 0;
r(1) +=
j > 0 ?
j * pow(x,
i) * pow(y,
j - 1) * pow(z,
k) : 0;
r(2) +=
k > 0 ?
k * pow(x,
i) * pow(y,
j) * pow(z,
k - 1) : 0;
}
}
}
}
}
};
const bool checkGradients;
:
DomainEleOp(
"FIELD1", OPROW), vAls(vals), diffVals(diff_vals),
checkGradients(check_grads) {}
const int nb_gauss_pts = getGaussPts().size2();
vAls.resize(nb_gauss_pts, false);
diffVals.resize(
SPACE_DIM, nb_gauss_pts,
false);
vAls.clear();
diffVals.clear();
}
if (nb_dofs) {
<<
"Type " << moab::CN::EntityTypeName(
type) <<
" side " << side;
for (int gg = 0; gg != nb_gauss_pts; gg++) {
for (int bb = 0; bb != nb_dofs; bb++) {
t_vals += t_base_fun * t_data;
++t_base_fun;
++t_data;
}
++t_vals;
}
if (checkGradients) {
auto t_diff_vals = getFTensor1FromMat<SPACE_DIM>(diffVals);
for (int gg = 0; gg != nb_gauss_pts; gg++) {
for (int bb = 0; bb != nb_dofs; bb++) {
t_diff_vals(
i) += t_diff_base_fun(
i) * t_data;
++t_diff_base_fun;
++t_data;
}
if (!std::isnormal(t_diff_vals(
d)))
++t_diff_vals;
}
}
}
}
};
boost::shared_ptr<VectorDouble> ptrVals;
boost::shared_ptr<MatrixDouble> ptrDiffVals;
const bool checkGradients;
boost::shared_ptr<VectorDouble> &ptr_vals,
boost::shared_ptr<MatrixDouble> &ptr_diff_vals,
bool check_grads)
:
DomainEleOp(
"FIELD1", OPROW), vAls(vals), diffVals(diff_vals),
ptrVals(ptr_vals), ptrDiffVals(ptr_diff_vals),
checkGradients(check_grads) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
const int nb_gauss_pts = data.
getN().size1();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
double err_val;
err_val = std::abs(t_vals - t_ptr_vals);
MOFEM_LOG(
"AT", Sev::noisy) <<
"Val op error " << err_val;
"Wrong value from operator %4.3e", err_val);
const double x = getCoordsAtGaussPts()(gg, 0);
const double y = getCoordsAtGaussPts()(gg, 1);
const double z = getCoordsAtGaussPts()(gg, 2);
err_val = std::fabs(delta_val * delta_val);
MOFEM_LOG(
"AT", Sev::verbose) << err_val <<
" : " << t_vals;
err_val);
++t_vals;
++t_ptr_vals;
}
if (checkGradients) {
auto t_diff_vals = getFTensor1FromMat<SPACE_DIM>(diffVals);
auto t_ptr_diff_vals = getFTensor1FromMat<SPACE_DIM>(*ptrDiffVals);
for (int gg = 0; gg != nb_gauss_pts; gg++) {
double err_diff_val;
t_delta_diff_val(
i) = t_diff_vals(
i) - t_ptr_diff_vals(
i);
err_diff_val = sqrt(t_delta_diff_val(
i) * t_delta_diff_val(
i));
MOFEM_LOG(
"AT", Sev::noisy) <<
"Diff val op error " << err_diff_val;
"Wrong derivatives from operator %4.3e", err_diff_val);
const double x = getCoordsAtGaussPts()(gg, 0);
const double y = getCoordsAtGaussPts()(gg, 1);
const double z = getCoordsAtGaussPts()(gg, 2);
t_delta_diff_val(
i) = t_diff_vals(
i) - t_diff_anal(
i);
err_diff_val = sqrt(t_delta_diff_val(
i) * t_delta_diff_val(
i));
<< "Diff val " << err_diff_val << " : "
<< sqrt(t_diff_vals(
i) * t_diff_vals(
i)) <<
" : "
<< t_diff_vals(0) << " (" << t_diff_anal(0) << ") "
<< t_diff_vals(1) << " (" << t_diff_anal(1) << ") "
<< t_diff_vals(2) << " (" << t_diff_anal(2) << ")";
else
<< "Diff val " << err_diff_val << " : "
<< sqrt(t_diff_vals(
i) * t_diff_vals(
i)) <<
" : "
<< t_diff_vals(0) << " (" << t_diff_anal(0) << ") "
<< t_diff_vals(1) << " (" << t_diff_anal(1) << ")";
<< getCoords()(3 * 1 + 0) - getCoords()(3 * 0 + 0);
<< getCoords()(3 * 1 + 1) - getCoords()(3 * 0 + 1);
<< getCoords()(3 * 1 + 2) - getCoords()(3 * 0 + 2);
MOFEM_LOG(
"AT", Sev::verbose) <<
"Diff val error " << err_diff_val;
"Wrong derivative of value %4.3e %4.3e", err_diff_val,
t_diff_anal.l2());
++t_diff_vals;
++t_ptr_diff_vals;
}
}
}
};
int main(
int argc,
char *argv[]) {
try {
DMType dm_name = "DMMOFEM";
auto core_log = logging::core::get();
core_log->add_sink(
simple->getAddBoundaryFE() =
true;
enum bases {
AINSWORTH,
AINSWORTH_LOBATTO,
DEMKOWICZ,
BERNSTEIN,
LASBASETOP
};
const char *list_bases[] = {"ainsworth", "ainsworth_lobatto", "demkowicz",
"bernstein"};
PetscBool flg;
PetscInt choice_base_value = AINSWORTH;
LASBASETOP, &choice_base_value, &flg);
if (flg != PETSC_TRUE)
if (choice_base_value == AINSWORTH)
if (choice_base_value == AINSWORTH_LOBATTO)
else if (choice_base_value == DEMKOWICZ)
else if (choice_base_value == BERNSTEIN)
enum spaces { H1SPACE, L2SPACE, LASBASETSPACE };
const char *list_spaces[] = {"h1", "l2"};
PetscInt choice_space_value = H1SPACE;
LASBASETSPACE, &choice_space_value, &flg);
if (flg != PETSC_TRUE)
if (choice_space_value == H1SPACE)
else if (choice_space_value == L2SPACE)
PETSC_NULL);
CHKERR moab.get_entities_by_dimension(0, 1, edges);
CHKERR moab.get_entities_by_dimension(0, 2, faces);
if (choice_base_value != BERNSTEIN) {
for (auto e : edges) {
rise_order.insert(e);
}
}
}
}
for (auto e : rise_order) {
rise_twice.insert(e);
}
}
}
std::vector<EntityHandle> &adjacency) {
switch (field.getSpace()) {
CHKERR moab.get_connectivity(&fe_ent, 1, adjacency,
true);
CHKERR moab.get_adjacencies(&fe_ent, 1, 1,
false, adjacency,
moab::Interface::UNION);
CHKERR moab.get_adjacencies(&fe_ent, 1, 2,
false, adjacency,
moab::Interface::UNION);
adjacency.push_back(fe_ent);
for (auto ent : adjacency)
fe.getSideNumberPtr(ent);
break;
CHKERR moab.get_entities_by_handle(field.getMeshset(), adjacency,
false);
for (auto ent : adjacency) {
.insert(
boost::shared_ptr<SideNumber>(
new SideNumber(ent, -1, 0, 0)));
}
} break;
default:
"this field is not implemented for TRI finite element");
}
};
simple->getDomainFEName(), MBTET, volume_adj);
simple->getDomainFEName(), MBHEX, volume_adj);
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
auto assemble_matrices_and_vectors = [&]() {
PETSC>::LinearForm<GAUSS>::OpSource<1, 1>;
pipeline_mng->getOpDomainLhsPipeline(), {NOSPACE});
pipeline_mng->getOpDomainRhsPipeline(), {NOSPACE});
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpMass(
"FIELD1", "FIELD1", [](double, double, double) { return 1.; }));
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
return 2 * p_data + 1;
};
};
auto solve_problem = [&] {
auto solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(solver);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
};
auto check_solution = [&] {
auto ptr_values = boost::make_shared<VectorDouble>();
auto ptr_diff_vals = boost::make_shared<MatrixDouble>();
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getOpDomainRhsPipeline().clear();
pipeline_mng->getOpDomainRhsPipeline(), {space});
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
ptr_diff_vals));
vals, diff_vals, ptr_values, ptr_diff_vals, true));
CHKERR pipeline_mng->loopFiniteElements();
};
auto post_proc = [&] {
auto get_domain_post_proc_fe = [&](auto post_proc_mesh_ptr) {
auto post_proc_fe =
boost::make_shared<PostProcBrokenMeshInMoabBaseCont<DomainEle>>(
m_field, post_proc_mesh_ptr);
post_proc_fe->getOpPtrVector(), {space});
auto ptr_values = boost::make_shared<VectorDouble>();
auto ptr_diff_vals = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
ptr_diff_vals));
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"FIELD1", ptr_values}},
{{"FIELD1_GRAD", ptr_diff_vals}},
{}, {})
);
return post_proc_fe;
};
auto get_bdy_post_proc_fe = [&](auto post_proc_mesh_ptr) {
auto bdy_post_proc_fe =
boost::make_shared<PostProcBrokenMeshInMoabBaseCont<BoundaryEle>>(
m_field, post_proc_mesh_ptr);
boost::make_shared<
auto ptr_values = boost::make_shared<VectorDouble>();
auto ptr_diff_vals = boost::make_shared<MatrixDouble>();
op_loop_side->getOpPtrVector(), {space});
op_loop_side->getOpPtrVector().push_back(
op_loop_side->getOpPtrVector().push_back(
ptr_diff_vals));
bdy_post_proc_fe->getOpPtrVector().push_back(op_loop_side);
bdy_post_proc_fe->getOpPtrVector().push_back(
bdy_post_proc_fe->getPostProcMesh(),
bdy_post_proc_fe->getMapGaussPts(),
{{"FIELD1", ptr_values}},
{{"FIELD1_GRAD", ptr_diff_vals}},
{}, {})
);
return bdy_post_proc_fe;
};
auto post_proc_mesh_ptr = boost::make_shared<moab::Core>();
auto post_proc_begin =
boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(
m_field, post_proc_mesh_ptr);
auto post_proc_end =
boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
m_field, post_proc_mesh_ptr);
auto domain_post_proc_fe = get_domain_post_proc_fe(post_proc_mesh_ptr);
auto bdy_post_proc_fe = get_bdy_post_proc_fe(post_proc_mesh_ptr);
domain_post_proc_fe);
bdy_post_proc_fe);
CHKERR post_proc_end->writeFile(
"out_post_proc.h5m");
};
CHKERR assemble_matrices_and_vectors();
}
}