284 {
285
287
288 try {
289
290 DMType dm_name = "DMMOFEM";
292
293 moab::Core mb_instance;
294 moab::Interface &moab = mb_instance;
295
296
297 auto core_log = logging::core::get();
298 core_log->add_sink(
302
303
306
310
311 simple->getAddBoundaryFE() =
true;
312
314
315
316 enum bases {
317 AINSWORTH,
318 AINSWORTH_LOBATTO,
319 DEMKOWICZ,
320 BERNSTEIN,
321 LASBASETOP
322 };
323 const char *list_bases[] = {"ainsworth", "ainsworth_labatto", "demkowicz",
324 "bernstein"};
325 PetscBool flg;
326 PetscInt choice_base_value = AINSWORTH;
328 LASBASETOP, &choice_base_value, &flg);
329
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)
341
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)
354
356 PETSC_NULL);
357
358 CHKERR simple->addDomainField(
"FIELD1", space, base, 1);
360
362 CHKERR moab.get_entities_by_dimension(0, 1, edges);
363 CHKERR moab.get_entities_by_dimension(0, 2, faces);
364
365 if (choice_base_value != BERNSTEIN) {
367
369 for (auto e : edges) {
371 rise_order.insert(e);
372 }
374 }
375
376 for (auto f : faces) {
378 rise_order.insert(f);
379 }
381 }
382
384 for (auto e : rise_order) {
386 rise_twice.insert(e);
387 }
389 }
390
392
394 }
395
397
398 auto volume_adj = [](moab::Interface &moab,
const Field &field,
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);
414
415 for (auto ent : adjacency)
416 fe.getSideNumberPtr(ent);
417 break;
419 CHKERR moab.get_entities_by_handle(field.getMeshset(), adjacency,
420 false);
421 for (auto ent : adjacency) {
423 .insert(
424 boost::shared_ptr<SideNumber>(
new SideNumber(ent, -1, 0, 0)));
425 }
426 } break;
427 default:
429 "this field is not implemented for TRI finite element");
430 }
431
433 };
434
436 simple->getDomainFEName(), MBTET, volume_adj);
438 simple->getDomainFEName(), MBHEX, volume_adj);
439
444
445 auto dm =
simple->getDM();
446
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>();
452
453 auto assemble_matrices_and_vectors = [&]() {
455
458
460 pipeline_mng->getOpDomainLhsPipeline(), {NOSPACE});
462 pipeline_mng->getOpDomainRhsPipeline(), {NOSPACE});
463
466
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(
473
474 pipeline_mng->getOpDomainRhsPipeline().push_back(
475 new OpSource("FIELD1", ApproxFunctions::fUn));
476
478 return 2 * p_data + 1;
479 };
482
484 };
485
486 auto solve_problem = [&] {
488 auto solver = pipeline_mng->createKSP();
489 CHKERR KSPSetFromOptions(solver);
491
492 auto dm =
simple->getDM();
495
497 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
498 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
500
502 };
503
504 auto check_solution = [&] {
506
507 auto ptr_values = boost::make_shared<VectorDouble>();
508 auto ptr_diff_vals = boost::make_shared<MatrixDouble>();
509
510 pipeline_mng->getDomainLhsFE().reset();
511 pipeline_mng->getOpDomainRhsPipeline().clear();
512
514 pipeline_mng->getOpDomainRhsPipeline(), {space});
515
516 pipeline_mng->getOpDomainRhsPipeline().push_back(
518 pipeline_mng->getOpDomainRhsPipeline().push_back(
520 pipeline_mng->getOpDomainRhsPipeline().push_back(
522 ptr_diff_vals));
524 vals, diff_vals, ptr_values, ptr_diff_vals, true));
525
526 CHKERR pipeline_mng->loopFiniteElements();
527
529 };
530
531 auto post_proc = [&] {
533
535
536 auto get_domain_post_proc_fe = [&](auto post_proc_mesh_ptr) {
537 auto post_proc_fe =
538 boost::make_shared<PostProcBrokenMeshInMoabBaseCont<DomainEle>>(
539 m_field, post_proc_mesh_ptr);
540
542 post_proc_fe->getOpPtrVector(), {space});
543
544 auto ptr_values = boost::make_shared<VectorDouble>();
545 auto ptr_diff_vals = boost::make_shared<MatrixDouble>();
546
547 post_proc_fe->getOpPtrVector().push_back(
549 post_proc_fe->getOpPtrVector().push_back(
551 ptr_diff_vals));
552
553 post_proc_fe->getOpPtrVector().push_back(
554
556
557 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
558
559 {{"FIELD1", ptr_values}},
560
561 {{"FIELD1_GRAD", ptr_diff_vals}},
562
563 {}, {})
564
565 );
566 return post_proc_fe;
567 };
568
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);
573
576
577 auto ptr_values = boost::make_shared<VectorDouble>();
578 auto ptr_diff_vals = boost::make_shared<MatrixDouble>();
579
580
582 op_loop_side->getOpPtrVector(), {space});
583 op_loop_side->getOpPtrVector().push_back(
585 op_loop_side->getOpPtrVector().push_back(
587 ptr_diff_vals));
588
589 bdy_post_proc_fe->getOpPtrVector().push_back(op_loop_side);
590
591 bdy_post_proc_fe->getOpPtrVector().push_back(
592
594
595 bdy_post_proc_fe->getPostProcMesh(),
596 bdy_post_proc_fe->getMapGaussPts(),
597
598 {{"FIELD1", ptr_values}},
599
600 {{"FIELD1_GRAD", ptr_diff_vals}},
601
602 {}, {})
603
604 );
605
606 return bdy_post_proc_fe;
607 };
608
609 auto post_proc_mesh_ptr = boost::make_shared<moab::Core>();
610 auto post_proc_begin =
611 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(
612 m_field, post_proc_mesh_ptr);
613 auto post_proc_end =
614 boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
615 m_field, post_proc_mesh_ptr);
616 auto domain_post_proc_fe = get_domain_post_proc_fe(post_proc_mesh_ptr);
617 auto bdy_post_proc_fe = get_bdy_post_proc_fe(post_proc_mesh_ptr);
618
621 domain_post_proc_fe);
623 bdy_post_proc_fe);
625 CHKERR post_proc_end->writeFile(
"out_post_proc.h5m");
626
628 };
629
630 CHKERR assemble_matrices_and_vectors();
634 }
636
638}
void simple(double P1[], double P2[], double P3[], double c[], const int N)
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
FieldSpace
approximation spaces
@ L2
field with C-1 continuity
@ NOFIELD
scalar or vector of scalars describe (no true field)
@ HCURL
field with continuous tangents
@ HDIV
field with continuous normal traction
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
multi_index_container< boost::shared_ptr< SideNumber >, indexed_by< ordered_unique< member< SideNumber, EntityHandle, &SideNumber::ent > >, ordered_non_unique< composite_key< SideNumber, const_mem_fun< SideNumber, EntityType, &SideNumber::getEntType >, member< SideNumber, signed char, &SideNumber::side_number > > > > > SideNumber_multiIndex
SideNumber_multiIndex for SideNumber.
virtual MoFEMErrorCode modify_finite_element_adjacency_table(const std::string &fe_name, const EntityType type, ElementAdjacencyFunct function)=0
modify finite element table, only for advanced user
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
FTensor::Index< 'i', SPACE_DIM > i
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Add operators pushing bases from local to physical configuration.
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Finite element data for entity.
Provide data structure for (tensor) field approximation.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Element used to execute operators on side of the element.
Post post-proc data at points from hash maps.
Clear Schur complement internal data.
Assemble Schur complement.
PipelineManager interface.
keeps information about side number for the finite element
Simple interface for fast problem set-up.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.