v0.14.0
helmholtz.cpp

Using PipelineManager interface calculate Helmholtz problem. Example show how to manage complex variable fields.

/**
* \file helmholtz.cpp
* \example helmholtz.cpp
*
* Using PipelineManager interface calculate Helmholtz problem. Example show how
* to manage complex variable fields.
*/
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
constexpr int SPACE_DIM = 2;
PETSC>::LinearForm<GAUSS>::OpSource<1, 1>;
struct Example {
Example(MoFEM::Interface &m_field) : mField(m_field) {}
MoFEMErrorCode runProblem();
private:
Simple *simpleInterface;
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
MoFEMErrorCode readMesh();
MoFEMErrorCode setupProblem();
MoFEMErrorCode boundaryCondition();
MoFEMErrorCode assembleSystem();
MoFEMErrorCode solveSystem();
MoFEMErrorCode outputResults();
MoFEMErrorCode checkResults();
};
//! [run problem]
CHKERR readMesh();
CHKERR setupProblem();
CHKERR boundaryCondition();
CHKERR assembleSystem();
CHKERR solveSystem();
CHKERR outputResults();
CHKERR checkResults();
}
//! [run problem]
//! [Read mesh]
CHKERR mField.getInterface(simpleInterface);
CHKERR simpleInterface->getOptions();
CHKERR simpleInterface->loadFile();
}
//! [Read mesh]
//! [Set up problem]
// Add field
CHKERR simpleInterface->addDomainField("P_REAL", H1,
CHKERR simpleInterface->addDomainField("P_IMAG", H1,
CHKERR simpleInterface->addBoundaryField("P_REAL", H1,
CHKERR simpleInterface->addBoundaryField("P_IMAG", H1,
int order = 6;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
CHKERR simpleInterface->setFieldOrder("P_REAL", order);
CHKERR simpleInterface->setFieldOrder("P_IMAG", order);
CHKERR simpleInterface->setUp();
}
//! [Set up problem]
//! [Applying essential BC]
auto get_ents_on_mesh_skin = [&]() {
Range boundary_entities;
std::string entity_name = it->getName();
if (entity_name.compare(0, 2, "BC") == 0) {
CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 1,
boundary_entities, true);
}
}
// Add vertices to boundary entities
Range boundary_vertices;
CHKERR mField.get_moab().get_connectivity(boundary_entities,
boundary_vertices, true);
boundary_entities.merge(boundary_vertices);
return boundary_entities;
};
auto mark_boundary_dofs = [&](Range &&skin_edges) {
auto problem_manager = mField.getInterface<ProblemsManager>();
auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
problem_manager->markDofs(simpleInterface->getProblemName(), ROW,
skin_edges, *marker_ptr);
return marker_ptr;
};
auto remove_dofs_from_problem = [&](Range &&skin_edges) {
auto problem_manager = mField.getInterface<ProblemsManager>();
CHKERR problem_manager->removeDofsOnEntities(
simpleInterface->getProblemName(), "P_IMAG", skin_edges, 0, 1);
};
// Get global local vector of marked DOFs. Is global, since is set for all
// DOFs on processor. Is local since only DOFs on processor are in the
// vector. To access DOFs use local indices.
boundaryMarker = mark_boundary_dofs(get_ents_on_mesh_skin());
CHKERR remove_dofs_from_problem(get_ents_on_mesh_skin());
}
//! [Applying essential BC]
//! [Push operators to pipeline]
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
double k = 90;
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-k", &k, PETSC_NULL);
auto beta = [](const double, const double, const double) { return -1; };
auto k2 = [k](const double, const double, const double) { return pow(k, 2); };
auto kp = [k](const double, const double, const double) { return k; };
auto km = [k](const double, const double, const double) { return -k; };
auto integration_rule = [](int, int, int p_data) { return 2 * p_data; };
auto set_domain = [&]() {
pipeline_mng->getOpDomainLhsPipeline(), {H1});
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpSetBc("P_REAL", true, boundaryMarker));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpDomainGradGrad("P_REAL", "P_REAL", beta));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpDomainGradGrad("P_IMAG", "P_IMAG", beta));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpDomainMass("P_REAL", "P_REAL", k2));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpDomainMass("P_IMAG", "P_IMAG", k2));
pipeline_mng->getOpDomainLhsPipeline().push_back(new OpUnSetBc("P_REAL"));
};
auto set_boundary = [&]() {
pipeline_mng->getOpBoundaryLhsPipeline(), {});
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
new OpSetBc("P_REAL", true, boundaryMarker));
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
new OpBoundaryMass("P_IMAG", "P_REAL", kp));
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
new OpBoundaryMass("P_REAL", "P_IMAG", km));
pipeline_mng->getOpBoundaryLhsPipeline().push_back(new OpUnSetBc("P_REAL"));
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
new OpSetBc("P_REAL", false, boundaryMarker));
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
new OpBoundaryMass("P_REAL", "P_REAL", beta));
pipeline_mng->getOpBoundaryLhsPipeline().push_back(new OpUnSetBc("P_REAL"));
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
new OpSetBc("P_REAL", false, boundaryMarker));
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
new OpBoundarySource("P_REAL", beta));
pipeline_mng->getOpBoundaryRhsPipeline().push_back(new OpUnSetBc("P_REAL"));
};
CHKERR set_domain();
CHKERR set_boundary();
}
//! [Push operators to pipeline]
//! [Solve]
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
auto solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(solver);
CHKERR KSPSetUp(solver);
auto dm = simpleInterface->getDM();
auto D = createDMVector(dm);
auto F = vectorDuplicate(D);
CHKERR KSPSolve(solver, F, D);
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
}
//! [Solve]
//! [Postprocess results]
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getDomainRhsFE().reset();
pipeline_mng->getBoundaryLhsFE().reset();
pipeline_mng->getBoundaryRhsFE().reset();
auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
auto p_real_ptr = boost::make_shared<VectorDouble>();
auto p_imag_ptr = boost::make_shared<VectorDouble>();
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("P_REAL", p_real_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("P_IMAG", p_imag_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"P_REAL", p_real_ptr}, {"P_IMAG", p_imag_ptr}},
{}, {}, {}
)
);
pipeline_mng->getDomainRhsFE() = post_proc_fe;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR post_proc_fe->writeFile("out_helmholtz.h5m");
}
//! [Postprocess results]
//! [Check results]
auto dm = simpleInterface->getDM();
auto D = createDMVector(dm);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_FORWARD);
double nrm2;
CHKERR VecNorm(D, NORM_2, &nrm2);
MOFEM_LOG("WORLD", Sev::inform)
<< std::setprecision(9) << "Solution norm " << nrm2;
PetscBool test_flg = PETSC_FALSE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test", &test_flg, PETSC_NULL);
if (test_flg) {
constexpr double regression_test = 97.261672;
constexpr double eps = 1e-6;
if (std::abs(nrm2 - regression_test) / regression_test > eps)
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Not converged solution");
}
}
//! [Check results]
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
MoFEM::Core::Initialize(&argc, &argv, param_file, help);
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
//! [Register MoFEM discrete manager in PETSc
//! [Create MoAB]
moab::Core mb_instance; ///< mesh database
moab::Interface &moab = mb_instance; ///< mesh database interface
//! [Create MoAB]
//! [Create MoFEM]
MoFEM::Core core(moab); ///< finite element database
MoFEM::Interface &m_field = core; ///< finite element database insterface
//! [Create MoFEM]
//! [Example]
Example ex(m_field);
CHKERR ex.runProblem();
//! [Example]
}
}
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::EdgeElementForcesAndSourcesCore
Edge finite element.
Definition: EdgeElementForcesAndSourcesCore.hpp:30
Example::checkResults
MoFEMErrorCode checkResults()
[Postprocess results]
Definition: dynamic_first_order_con_law.cpp:1205
SPACE_DIM
constexpr int SPACE_DIM
Definition: helmholtz.cpp:22
main
int main(int argc, char *argv[])
[Check results]
Definition: helmholtz.cpp:313
MoFEM::PipelineManager::getDomainRhsFE
boost::shared_ptr< FEMethod > & getDomainRhsFE()
Definition: PipelineManager.hpp:405
Example::assembleSystem
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
Definition: dynamic_first_order_con_law.cpp:647
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
MoFEM::PipelineManager::loopFiniteElements
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
Definition: PipelineManager.cpp:19
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::PipelineManager::getBoundaryLhsFE
boost::shared_ptr< FEMethod > & getBoundaryLhsFE()
Definition: PipelineManager.hpp:409
Example::simpleInterface
Simple * simpleInterface
Definition: helmholtz.cpp:41
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM.hpp
MoFEM::DMoFEMMeshToLocalVector
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:523
MoFEM::PipelineManager::getBoundaryRhsFE
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
Definition: PipelineManager.hpp:413
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
OpDomainMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
Definition: child_and_parent.cpp:53
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
OpBoundaryMass
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: seepage.cpp:71
ROW
@ ROW
Definition: definitions.h:136
MoFEM::PipelineManager::getOpBoundaryLhsPipeline
boost::ptr_deque< UserDataOperator > & getOpBoundaryLhsPipeline()
Get the Op Boundary Lhs Pipeline object.
Definition: PipelineManager.hpp:797
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
Example::readMesh
MoFEMErrorCode readMesh()
[Run problem]
Definition: dynamic_first_order_con_law.cpp:463
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:667
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:747
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
Example
[Example]
Definition: plastic.cpp:177
double
OpBoundarySource
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
Definition: helmholtz.cpp:31
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:317
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
Example::solveSystem
MoFEMErrorCode solveSystem()
[Solve]
Definition: dynamic_first_order_con_law.cpp:893
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
help
static char help[]
Definition: helmholtz.cpp:13
MoFEM::PipelineManager::createKSP
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
Definition: PipelineManager.cpp:49
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::ProblemsManager::markDofs
MoFEMErrorCode markDofs(const std::string problem_name, RowColData rc, const enum MarkOP op, const Range ents, std::vector< unsigned char > &marker) const
Create vector with marked indices.
Definition: ProblemsManager.cpp:3548
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
MoFEM::PipelineManager::setBoundaryLhsIntegrationRule
MoFEMErrorCode setBoundaryLhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:557
MoFEM::PipelineManager::getOpDomainLhsPipeline
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
Definition: PipelineManager.hpp:749
MoFEM::ForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: ForcesAndSourcesCore.hpp:482
BiLinearForm
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:413
MoFEM::EdgeElementForcesAndSourcesCore::UserDataOperator
default operator for EDGE element
Definition: EdgeElementForcesAndSourcesCore.hpp:68
MoFEM::PipelineManager::setDomainLhsIntegrationRule
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:503
integration_rule
auto integration_rule
Definition: free_surface.cpp:185
MoFEM::OpUnSetBc
Definition: FormsIntegrators.hpp:49
Range
MoFEM::PipelineManager::getDomainLhsFE
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Definition: PipelineManager.hpp:401
DomainEleOp
MoFEM::CoreTmp< 0 >::Initialize
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
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
OpDomainGradGrad
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
Definition: helmholtz.cpp:25
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
Definition: MeshsetsManager.hpp:71
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
Example::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:225
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
Example::boundaryCondition
MoFEMErrorCode boundaryCondition()
[Set up problem]
Definition: dynamic_first_order_con_law.cpp:536
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
MoFEM::OpSetBc
Set indices on entities on finite element.
Definition: FormsIntegrators.hpp:38
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::PipelineManager::getOpBoundaryRhsPipeline
boost::ptr_deque< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
Definition: PipelineManager.hpp:821
Example::runProblem
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:213
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
MoFEM::PipelineManager::setBoundaryRhsIntegrationRule
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:584
Example::mField
MoFEM::Interface & mField
Definition: plastic.cpp:184
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
Example::outputResults
MoFEMErrorCode outputResults()
[Solve]
Definition: dynamic_first_order_con_law.cpp:1183
convert.int
int
Definition: convert.py:64
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
OpDomainMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpDomainMass
Definition: helmholtz.cpp:27
OpBoundaryMass
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpBoundaryMass
Definition: helmholtz.cpp:29
F
@ F
Definition: free_surface.cpp:394
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698