v0.14.0
Loading...
Searching...
No Matches
Classes | Typedefs | Functions | Variables
scalar_check_approximation.cpp File Reference
#include <MoFEM.hpp>

Go to the source code of this file.

Classes

struct  ApproxFunctionsImpl< DIM >
 
struct  ElementsAndOps< DIM >
 
struct  ApproxFunctionsImpl< 2 >
 
struct  ApproxFunctionsImpl< 3 >
 
struct  OpValsDiffVals
 
struct  OpCheckValsDiffVals
 

Typedefs

using DomainEle = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::DomainEle
 
using BoundaryEle
 
using EleOnSide = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle
 
using EntData = EntitiesFieldData::EntData
 
using DomainEleOp = DomainEle::UserDataOperator
 
using ApproxFunctions = ApproxFunctionsImpl<SPACE_DIM>
 

Functions

int main (int argc, char *argv[])
 

Variables

static char help [] = "...\n\n"
 
static int approx_order = 4
 
constexpr int SPACE_DIM
 

Typedef Documentation

◆ ApproxFunctions

using ApproxFunctions = ApproxFunctionsImpl<SPACE_DIM>

Definition at line 96 of file scalar_check_approximation.cpp.

◆ BoundaryEle

using BoundaryEle

◆ DomainEle

Definition at line 24 of file scalar_check_approximation.cpp.

◆ DomainEleOp

Definition at line 30 of file scalar_check_approximation.cpp.

◆ EleOnSide

Definition at line 27 of file scalar_check_approximation.cpp.

◆ EntData

Definition at line 29 of file scalar_check_approximation.cpp.

Function Documentation

◆ main()

int main ( int argc,
char * argv[] )

Definition at line 284 of file scalar_check_approximation.cpp.

284 {
285
286 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
287
288 try {
289
290 DMType dm_name = "DMMOFEM";
291 CHKERR DMRegister_MoFEM(dm_name);
292
293 moab::Core mb_instance;
294 moab::Interface &moab = mb_instance;
295
296 // Add logging channel for example
297 auto core_log = logging::core::get();
298 core_log->add_sink(
300 LogManager::setLog("AT");
301 MOFEM_LOG_TAG("AT", "atom_test");
302
303 // Create MoFEM instance
304 MoFEM::Core core(moab);
305 MoFEM::Interface &m_field = core;
306
307 Simple *simple = m_field.getInterface<Simple>();
308 PipelineManager *pipeline_mng = m_field.getInterface<PipelineManager>();
309 CHKERR simple->getOptions();
310
311 simple->getAddBoundaryFE() = true;
312
313 CHKERR simple->loadFile("", "");
314
315 // Declare elements
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;
327 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
328 LASBASETOP, &choice_base_value, &flg);
329
330 if (flg != PETSC_TRUE)
331 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "base not set");
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;
345 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-space", list_spaces,
346 LASBASETSPACE, &choice_space_value, &flg);
347 if (flg != PETSC_TRUE)
348 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "space not set");
349 FieldSpace space = H1;
350 if (choice_space_value == H1SPACE)
351 space = H1;
352 else if (choice_space_value == L2SPACE)
353 space = L2;
354
355 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &approx_order,
356 PETSC_NULL);
357
358 CHKERR simple->addDomainField("FIELD1", space, base, 1);
359 CHKERR simple->setFieldOrder("FIELD1", approx_order);
360
361 Range edges, faces;
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) {
366 Range rise_order;
367
368 int i = 0;
369 for (auto e : edges) {
370 if (!(i % 2)) {
371 rise_order.insert(e);
372 }
373 ++i;
374 }
375
376 for (auto f : faces) {
377 if (!(i % 3)) {
378 rise_order.insert(f);
379 }
380 ++i;
381 }
382
383 Range rise_twice;
384 for (auto e : rise_order) {
385 if (!(i % 2)) {
386 rise_twice.insert(e);
387 }
388 ++i;
389 }
390
391 CHKERR simple->setFieldOrder("FIELD1", approx_order + 1, &rise_order);
392
393 CHKERR simple->setFieldOrder("FIELD1", approx_order + 2, &rise_twice);
394 }
395
396 CHKERR simple->defineFiniteElements();
397
398 auto volume_adj = [](moab::Interface &moab, const Field &field,
399 const EntFiniteElement &fe,
400 std::vector<EntityHandle> &adjacency) {
402 EntityHandle fe_ent = fe.getEnt();
403 switch (field.getSpace()) {
404 case H1:
405 CHKERR moab.get_connectivity(&fe_ent, 1, adjacency, true);
406 case HCURL:
407 CHKERR moab.get_adjacencies(&fe_ent, 1, 1, false, adjacency,
408 moab::Interface::UNION);
409 case HDIV:
410 case L2:
411 CHKERR moab.get_adjacencies(&fe_ent, 1, 2, false, adjacency,
412 moab::Interface::UNION);
413 adjacency.push_back(fe_ent);
414 // build side table
415 for (auto ent : adjacency)
416 fe.getSideNumberPtr(ent);
417 break;
418 case NOFIELD: {
419 CHKERR moab.get_entities_by_handle(field.getMeshset(), adjacency,
420 false);
421 for (auto ent : adjacency) {
422 const_cast<SideNumber_multiIndex &>(fe.getSideNumberTable())
423 .insert(
424 boost::shared_ptr<SideNumber>(new SideNumber(ent, -1, 0, 0)));
425 }
426 } break;
427 default:
428 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
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
440 CHKERR simple->defineProblem(PETSC_TRUE);
441 CHKERR simple->buildFields();
442 CHKERR simple->buildFiniteElements();
443 CHKERR simple->buildProblem();
444
445 auto dm = simple->getDM();
446
447 VectorDouble vals;
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>();
451 MatrixDouble diff_vals;
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(
472 new OpSchurAssembleEnd<SCHUR_DSYSV>({}, {}, {}, {}, {}));
473
474 pipeline_mng->getOpDomainRhsPipeline().push_back(
475 new OpSource("FIELD1", ApproxFunctions::fUn));
476
477 auto integration_rule = [](int, int, int p_data) {
478 return 2 * p_data + 1;
479 };
480 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
481 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
482
484 };
485
486 auto solve_problem = [&] {
488 auto solver = pipeline_mng->createKSP();
489 CHKERR KSPSetFromOptions(solver);
490 CHKERR KSPSetUp(solver);
491
492 auto dm = simple->getDM();
493 auto D = createDMVector(dm);
494 auto F = vectorDuplicate(D);
495
496 CHKERR KSPSolve(solver, F, D);
497 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
498 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
499 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
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(
517 new OpValsDiffVals(vals, diff_vals, true));
518 pipeline_mng->getOpDomainRhsPipeline().push_back(
519 new OpCalculateScalarFieldValues("FIELD1", ptr_values));
520 pipeline_mng->getOpDomainRhsPipeline().push_back(
522 ptr_diff_vals));
523 pipeline_mng->getOpDomainRhsPipeline().push_back(new OpCheckValsDiffVals(
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(
548 new OpCalculateScalarFieldValues("FIELD1", ptr_values));
549 post_proc_fe->getOpPtrVector().push_back(
551 ptr_diff_vals));
552
553 post_proc_fe->getOpPtrVector().push_back(
554
555 new OpPPMap(
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
574 auto op_loop_side = new OpLoopSide<EleOnSide>(
575 m_field, simple->getDomainFEName(), SPACE_DIM, Sev::noisy,
576 boost::make_shared<
577 ForcesAndSourcesCore::UserDataOperator::AdjCache>());
578
579 auto ptr_values = boost::make_shared<VectorDouble>();
580 auto ptr_diff_vals = boost::make_shared<MatrixDouble>();
581
582 // push operators to side element
584 op_loop_side->getOpPtrVector(), {space});
585 op_loop_side->getOpPtrVector().push_back(
586 new OpCalculateScalarFieldValues("FIELD1", ptr_values));
587 op_loop_side->getOpPtrVector().push_back(
589 ptr_diff_vals));
590 // push op to boundary element
591 bdy_post_proc_fe->getOpPtrVector().push_back(op_loop_side);
592
593 bdy_post_proc_fe->getOpPtrVector().push_back(
594
595 new OpPPMap(
596
597 bdy_post_proc_fe->getPostProcMesh(),
598 bdy_post_proc_fe->getMapGaussPts(),
599
600 {{"FIELD1", ptr_values}},
601
602 {{"FIELD1_GRAD", ptr_diff_vals}},
603
604 {}, {})
605
606 );
607
608 return bdy_post_proc_fe;
609 };
610
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);
615 auto post_proc_end =
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);
620
621 CHKERR DMoFEMPreProcessFiniteElements(dm, post_proc_begin->getFEMethod());
622 CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
623 domain_post_proc_fe);
624 CHKERR DMoFEMLoopFiniteElements(dm, simple->getBoundaryFEName(),
625 bdy_post_proc_fe);
626 CHKERR DMoFEMPostProcessFiniteElements(dm, post_proc_end->getFEMethod());
627 CHKERR post_proc_end->writeFile("out_post_proc.h5m");
628
630 };
631
632 CHKERR assemble_matrices_and_vectors();
633 CHKERR solve_problem();
634 CHKERR check_solution();
635 CHKERR post_proc();
636 }
638
640}
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
Definition definitions.h:58
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
@ AINSWORTH_LOBATTO_BASE
Definition definitions.h:62
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition definitions.h:64
FieldSpace
approximation spaces
Definition definitions.h:82
@ L2
field with C-1 continuity
Definition definitions.h:88
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition definitions.h:84
@ H1
continuous field
Definition definitions.h:85
@ HCURL
field with continuous tangents
Definition definitions.h:86
@ HDIV
field with continuous normal traction
Definition definitions.h:87
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_IMPOSSIBLE_CASE
Definition definitions.h:35
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
#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
@ F
auto integration_rule
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:542
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition DMMoFEM.cpp:509
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:47
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:572
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1003
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:532
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
double D
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(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
static int approx_order
static char help[]
constexpr int SPACE_DIM
Add operators pushing bases from local to physical configuration.
Core (interface) class.
Definition Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:112
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.
Definition Schur.hpp:22
Assemble Schur complement.
Definition Schur.hpp:104
PipelineManager interface.
keeps information about side number for the finite element
Simple interface for fast problem set-up.
Definition Simple.hpp:27
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.

Variable Documentation

◆ approx_order

int approx_order = 4
static

◆ help

char help[] = "...\n\n"
static

Definition at line 13 of file scalar_check_approximation.cpp.

◆ SPACE_DIM

constexpr int SPACE_DIM
constexpr
Initial value:
=
#define EXECUTABLE_DIMENSION
Definition plastic.cpp:13
Examples
scalar_check_approximation.cpp.

Definition at line 21 of file scalar_check_approximation.cpp.