11 using namespace MoFEM;
13 static char help[] =
"...\n\n";
33 static double fUn(
const double x,
const double y,
double z) {
36 for (
int i = 0;
i <= o; ++
i) {
39 r += pow(x,
i) * pow(y,
j);
49 for (
int i = 0;
i <= o; ++
i) {
52 r(0) +=
i > 0 ?
i * pow(x,
i - 1) * pow(y,
j) : 0;
53 r(1) +=
j > 0 ?
j * pow(x,
i) * pow(y,
j - 1) : 0;
62 static double fUn(
const double x,
const double y,
double z) {
65 for (
int i = 0;
i <= o; ++
i) {
66 for (
int j = 0;
j <= o -
i;
j++) {
69 r += pow(x,
i) * pow(y,
j) * pow(z,
k);
81 for (
int i = 0;
i <= o; ++
i) {
82 for (
int j = 0;
j <= o -
i;
j++) {
85 r(0) +=
i > 0 ?
i * pow(x,
i - 1) * pow(y,
j) * pow(z,
k) : 0;
86 r(1) +=
j > 0 ?
j * pow(x,
i) * pow(y,
j - 1) * pow(z,
k) : 0;
87 r(2) +=
k > 0 ?
k * pow(x,
i) * pow(y,
j) * pow(z,
k - 1) : 0;
103 :
DomainEleOp(
"FIELD1", OPROW), vAls(vals), diffVals(diff_vals),
104 checkGradients(check_grads) {}
111 const int nb_gauss_pts = getGaussPts().size2();
112 if (
type == MBVERTEX) {
113 vAls.resize(nb_gauss_pts,
false);
114 diffVals.resize(
SPACE_DIM, nb_gauss_pts,
false);
123 <<
"Type " << moab::CN::EntityTypeName(
type) <<
" side " << side;
129 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
131 for (
int bb = 0; bb != nb_dofs; bb++) {
132 t_vals += t_base_fun * t_data;
136 const double v = t_vals;
137 if (!std::isnormal(
v))
142 if (checkGradients) {
143 auto t_diff_vals = getFTensor1FromMat<SPACE_DIM>(diffVals);
145 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
147 for (
int bb = 0; bb != nb_dofs; bb++) {
148 t_diff_vals(
i) += t_diff_base_fun(
i) * t_data;
153 if (!std::isnormal(t_diff_vals(
d)))
171 boost::shared_ptr<VectorDouble> &ptr_vals,
172 boost::shared_ptr<MatrixDouble> &ptr_diff_vals,
174 :
DomainEleOp(
"FIELD1", OPROW), vAls(vals), diffVals(diff_vals),
175 ptrVals(ptr_vals), ptrDiffVals(ptr_diff_vals),
176 checkGradients(check_grads) {
177 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE],
false);
185 const double eps = 1e-6;
186 const int nb_gauss_pts = data.
getN().size1();
192 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
197 err_val = std::abs(t_vals - t_ptr_vals);
198 MOFEM_LOG(
"AT", Sev::noisy) <<
"Val op error " << err_val;
202 "Wrong value from operator %4.3e", err_val);
204 const double x = getCoordsAtGaussPts()(gg, 0);
205 const double y = getCoordsAtGaussPts()(gg, 1);
206 const double z = getCoordsAtGaussPts()(gg, 2);
211 err_val = std::fabs(delta_val * delta_val);
212 MOFEM_LOG(
"AT", Sev::verbose) << err_val <<
" : " << t_vals;
221 if (checkGradients) {
223 auto t_diff_vals = getFTensor1FromMat<SPACE_DIM>(diffVals);
224 auto t_ptr_diff_vals = getFTensor1FromMat<SPACE_DIM>(*ptrDiffVals);
226 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
231 t_delta_diff_val(
i) = t_diff_vals(
i) - t_ptr_diff_vals(
i);
232 err_diff_val = sqrt(t_delta_diff_val(
i) * t_delta_diff_val(
i));
233 MOFEM_LOG(
"AT", Sev::noisy) <<
"Diff val op error " << err_diff_val;
235 if (err_diff_val >
eps)
237 "Wrong derivatives from operator %4.3e", err_diff_val);
239 const double x = getCoordsAtGaussPts()(gg, 0);
240 const double y = getCoordsAtGaussPts()(gg, 1);
241 const double z = getCoordsAtGaussPts()(gg, 2);
245 t_delta_diff_val(
i) = t_diff_vals(
i) - t_diff_anal(
i);
247 err_diff_val = sqrt(t_delta_diff_val(
i) * t_delta_diff_val(
i));
250 <<
"Diff val " << err_diff_val <<
" : "
251 << sqrt(t_diff_vals(
i) * t_diff_vals(
i)) <<
" : "
252 << t_diff_vals(0) <<
" (" << t_diff_anal(0) <<
") "
253 << t_diff_vals(1) <<
" (" << t_diff_anal(1) <<
") "
254 << t_diff_vals(2) <<
" (" << t_diff_anal(2) <<
")";
257 <<
"Diff val " << err_diff_val <<
" : "
258 << sqrt(t_diff_vals(
i) * t_diff_vals(
i)) <<
" : "
259 << t_diff_vals(0) <<
" (" << t_diff_anal(0) <<
") "
260 << t_diff_vals(1) <<
" (" << t_diff_anal(1) <<
")";
263 << getCoords()(3 * 1 + 0) - getCoords()(3 * 0 + 0);
265 << getCoords()(3 * 1 + 1) - getCoords()(3 * 0 + 1);
267 << getCoords()(3 * 1 + 2) - getCoords()(3 * 0 + 2);
269 MOFEM_LOG(
"AT", Sev::verbose) <<
"Diff val error " << err_diff_val;
270 if (err_diff_val >
eps)
272 "Wrong derivative of value %4.3e %4.3e", err_diff_val,
284 int main(
int argc,
char *argv[]) {
290 DMType dm_name =
"DMMOFEM";
297 auto core_log = logging::core::get();
311 simple->getAddBoundaryFE() =
true;
323 const char *list_bases[] = {
"ainsworth",
"ainsworth_labatto",
"demkowicz",
326 PetscInt choice_base_value = AINSWORTH;
328 LASBASETOP, &choice_base_value, &flg);
330 if (flg != PETSC_TRUE)
333 if (choice_base_value == AINSWORTH)
335 if (choice_base_value == AINSWORTH_LOBATTO)
337 else if (choice_base_value == DEMKOWICZ)
339 else if (choice_base_value == BERNSTEIN)
342 enum spaces { H1SPACE, L2SPACE, LASBASETSPACE };
343 const char *list_spaces[] = {
"h1",
"l2"};
344 PetscInt choice_space_value = H1SPACE;
346 LASBASETSPACE, &choice_space_value, &flg);
347 if (flg != PETSC_TRUE)
350 if (choice_space_value == H1SPACE)
352 else if (choice_space_value == L2SPACE)
358 CHKERR simple->addDomainField(
"FIELD1", space, base, 1);
362 CHKERR moab.get_entities_by_dimension(0, 1, edges);
363 CHKERR moab.get_entities_by_dimension(0, 2, faces);
365 if (choice_base_value != BERNSTEIN) {
369 for (
auto e : edges) {
371 rise_order.insert(e);
376 for (
auto f : faces) {
378 rise_order.insert(
f);
384 for (
auto e : rise_order) {
386 rise_twice.insert(e);
400 std::vector<EntityHandle> &adjacency) {
403 switch (field.getSpace()) {
405 CHKERR moab.get_connectivity(&fe_ent, 1, adjacency,
true);
407 CHKERR moab.get_adjacencies(&fe_ent, 1, 1,
false, adjacency,
408 moab::Interface::UNION);
411 CHKERR moab.get_adjacencies(&fe_ent, 1, 2,
false, adjacency,
412 moab::Interface::UNION);
413 adjacency.push_back(fe_ent);
415 for (
auto ent : adjacency)
416 fe.getSideNumberPtr(ent);
419 CHKERR moab.get_entities_by_handle(field.getMeshset(), adjacency,
421 for (
auto ent : adjacency) {
424 boost::shared_ptr<SideNumber>(
new SideNumber(ent, -1, 0, 0)));
429 "this field is not implemented for TRI finite element");
436 simple->getDomainFEName(), MBTET, volume_adj);
438 simple->getDomainFEName(), MBHEX, volume_adj);
445 auto dm =
simple->getDM();
448 auto jac_ptr = boost::make_shared<MatrixDouble>();
449 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
450 auto det_ptr = boost::make_shared<VectorDouble>();
453 auto assemble_matrices_and_vectors = [&]() {
457 PETSC>::LinearForm<GAUSS>::OpSource<1, 1>;
460 pipeline_mng->getOpDomainLhsPipeline(), {NOSPACE});
462 pipeline_mng->getOpDomainRhsPipeline(), {NOSPACE});
467 pipeline_mng->getOpDomainLhsPipeline().push_back(
469 pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpMass(
470 "FIELD1",
"FIELD1", [](
double,
double,
double) {
return 1.; }));
471 pipeline_mng->getOpDomainLhsPipeline().push_back(
474 pipeline_mng->getOpDomainRhsPipeline().push_back(
478 return 2 * p_data + 1;
486 auto solve_problem = [&] {
488 auto solver = pipeline_mng->createKSP();
489 CHKERR KSPSetFromOptions(solver);
492 auto dm =
simple->getDM();
497 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
498 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
504 auto check_solution = [&] {
507 auto ptr_values = boost::make_shared<VectorDouble>();
508 auto ptr_diff_vals = boost::make_shared<MatrixDouble>();
510 pipeline_mng->getDomainLhsFE().reset();
511 pipeline_mng->getOpDomainRhsPipeline().clear();
514 pipeline_mng->getOpDomainRhsPipeline(), {space});
516 pipeline_mng->getOpDomainRhsPipeline().push_back(
518 pipeline_mng->getOpDomainRhsPipeline().push_back(
520 pipeline_mng->getOpDomainRhsPipeline().push_back(
524 vals, diff_vals, ptr_values, ptr_diff_vals,
true));
526 CHKERR pipeline_mng->loopFiniteElements();
531 auto post_proc = [&] {
536 auto get_domain_post_proc_fe = [&](
auto post_proc_mesh_ptr) {
538 boost::make_shared<PostProcBrokenMeshInMoabBaseCont<DomainEle>>(
539 m_field, post_proc_mesh_ptr);
542 post_proc_fe->getOpPtrVector(), {space});
544 auto ptr_values = boost::make_shared<VectorDouble>();
545 auto ptr_diff_vals = boost::make_shared<MatrixDouble>();
547 post_proc_fe->getOpPtrVector().push_back(
549 post_proc_fe->getOpPtrVector().push_back(
553 post_proc_fe->getOpPtrVector().push_back(
557 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
559 {{
"FIELD1", ptr_values}},
561 {{
"FIELD1_GRAD", ptr_diff_vals}},
569 auto get_bdy_post_proc_fe = [&](
auto post_proc_mesh_ptr) {
570 auto bdy_post_proc_fe =
571 boost::make_shared<PostProcBrokenMeshInMoabBaseCont<BoundaryEle>>(
572 m_field, post_proc_mesh_ptr);
579 auto ptr_values = boost::make_shared<VectorDouble>();
580 auto ptr_diff_vals = boost::make_shared<MatrixDouble>();
584 op_loop_side->getOpPtrVector(), {space});
585 op_loop_side->getOpPtrVector().push_back(
587 op_loop_side->getOpPtrVector().push_back(
591 bdy_post_proc_fe->getOpPtrVector().push_back(op_loop_side);
593 bdy_post_proc_fe->getOpPtrVector().push_back(
597 bdy_post_proc_fe->getPostProcMesh(),
598 bdy_post_proc_fe->getMapGaussPts(),
600 {{
"FIELD1", ptr_values}},
602 {{
"FIELD1_GRAD", ptr_diff_vals}},
608 return bdy_post_proc_fe;
611 auto post_proc_mesh_ptr = boost::make_shared<moab::Core>();
612 auto post_proc_begin =
613 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(
614 m_field, post_proc_mesh_ptr);
616 boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
617 m_field, post_proc_mesh_ptr);
618 auto domain_post_proc_fe = get_domain_post_proc_fe(post_proc_mesh_ptr);
619 auto bdy_post_proc_fe = get_bdy_post_proc_fe(post_proc_mesh_ptr);
623 domain_post_proc_fe);
627 CHKERR post_proc_end->writeFile(
"out_post_proc.h5m");
632 CHKERR assemble_matrices_and_vectors();