v0.14.0
Loading...
Searching...
No Matches
poisson_2d_homogeneous.cpp

Solution of poisson equation. Direct implementation of User Data Operators for teaching proposes.

Note
In practical application we suggest use form integrators to generalise and simplify code. However, here we like to expose user to ways how to implement data operator from scratch.
/**
* \file poisson_2d_homogeneous.cpp
* \example poisson_2d_homogeneous.cpp
*
* Solution of poisson equation. Direct implementation of User Data Operators
* for teaching proposes.
*
* \note In practical application we suggest use form integrators to generalise
* and simplify code. However, here we like to expose user to ways how to
* implement data operator from scratch.
*/
constexpr auto field_name = "U";
constexpr int SPACE_DIM =
EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
using namespace MoFEM;
static char help[] = "...\n\n";
public:
// Declaration of the main function to run analysis
private:
// Declaration of other main functions called in runProgram()
// MoFEM interfaces
// Field name and approximation order
int oRder;
};
: mField(m_field) {}
//! [Read mesh]
CHKERR simpleInterface->getOptions();
CHKERR simpleInterface->loadFile();
}
//! [Read mesh]
//! [Setup problem]
CHKERR simpleInterface->addDomainField(field_name, H1,
int oRder = 3;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &oRder, PETSC_NULL);
}
//! [Setup problem]
//! [Boundary condition]
auto bc_mng = mField.getInterface<BcManager>();
// Remove BCs from blockset name "BOUNDARY_CONDITION" or SETU, note that you
// can use regular expression to put list of blocksets;
CHKERR bc_mng->removeBlockDOFsOnEntities<BcScalarMeshsetType<BLOCKSET>>(
simpleInterface->getProblemName(), "(BOUNDARY_CONDITION|SETU)",
std::string(field_name), true);
}
//! [Boundary condition]
//! [Assemble system]
auto pipeline_mng = mField.getInterface<PipelineManager>();
{ // Push operators to the Pipeline that is responsible for calculating LHS
pipeline_mng->getOpDomainLhsPipeline(), {H1});
pipeline_mng->getOpDomainLhsPipeline().push_back(
}
{ // Push operators to the Pipeline that is responsible for calculating RHS
auto set_values_to_bc_dofs = [&](auto &fe) {
auto get_bc_hook = [&]() {
return hook;
};
fe->preProcessHook = get_bc_hook();
};
// you can skip that if boundary condition is prescribing zero
auto calculate_residual_from_set_values_on_bc = [&](auto &pipeline) {
using DomainEle =
auto grad_u_vals_ptr = boost::make_shared<MatrixDouble>();
pipeline_mng->getOpDomainRhsPipeline().push_back(
grad_u_vals_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpInternal(field_name, grad_u_vals_ptr,
[](double, double, double) constexpr { return -1; }));
};
pipeline_mng->getOpDomainRhsPipeline(), {H1});
set_values_to_bc_dofs(pipeline_mng->getDomainRhsFE());
calculate_residual_from_set_values_on_bc(
pipeline_mng->getOpDomainRhsPipeline());
pipeline_mng->getOpDomainRhsPipeline().push_back(
}
}
//! [Assemble system]
//! [Set integration rules]
auto rule_lhs = [](int, int, int p) -> int { return 2 * (p - 1); };
auto rule_rhs = [](int, int, int p) -> int { return p; };
auto pipeline_mng = mField.getInterface<PipelineManager>();
CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule_lhs);
CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule_rhs);
}
//! [Set integration rules]
//! [Solve system]
auto pipeline_mng = mField.getInterface<PipelineManager>();
auto ksp_solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(ksp_solver);
CHKERR KSPSetUp(ksp_solver);
// Create RHS and solution vectors
auto dm = simpleInterface->getDM();
auto F = createDMVector(dm);
auto D = vectorDuplicate(F);
// Solve the system
CHKERR KSPSolve(ksp_solver, F, D);
// Scatter result data on the mesh
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
}
//! [Solve system]
//! [Output results]
auto pipeline_mng = mField.getInterface<PipelineManager>();
pipeline_mng->getDomainLhsFE().reset();
auto post_proc_fe = boost::make_shared<PostProcFaceEle>(mField);
post_proc_fe->getOpPtrVector(), {H1});
auto u_ptr = boost::make_shared<VectorDouble>();
auto grad_u_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
OpPPMap::DataMapVec{{"U", u_ptr}},
OpPPMap::DataMapMat{{"GRAD_U", grad_u_ptr}},
OpPPMap::DataMapMat{},
OpPPMap::DataMapMat{}
)
);
pipeline_mng->getDomainRhsFE() = post_proc_fe;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR post_proc_fe->writeFile("out_result.h5m");
}
//! [Output results]
//! [Run program]
}
//! [Run program]
//! [Main]
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
// Error handling
try {
// Register MoFEM discrete manager in PETSc
DMType dm_name = "DMMOFEM";
// Create MOAB instance
moab::Core mb_instance; // mesh database
moab::Interface &moab = mb_instance; // mesh database interface
// Create MoFEM instance
MoFEM::Core core(moab); // finite element database
MoFEM::Interface &m_field = core; // finite element interface
// Run the main analysis
Poisson2DHomogeneous poisson_problem(m_field);
CHKERR poisson_problem.runProgram();
}
// Finish work: cleaning memory, getting statistics, etc.
return 0;
}
//! [Main]
static Index< 'p', 3 > p
std::string param_file
static char help[]
int main()
Definition: adol-c_atom.cpp:46
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
@ F
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
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1003
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
double D
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
constexpr auto field_name
constexpr int SPACE_DIM
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
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.
Class (Function) to enforce essential constrains.
Definition: Essential.hpp:25
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode boundaryCondition()
[Setup problem]
Poisson2DHomogeneous(MoFEM::Interface &m_field)
MoFEMErrorCode runProgram()
[Output results]
MoFEMErrorCode readMesh()
[Read mesh]
MoFEMErrorCode outputResults()
[Solve system]
MoFEMErrorCode assembleSystem()
[Boundary condition]
MoFEMErrorCode solveSystem()
[Set integration rules]
MoFEMErrorCode setIntegrationRules()
[Assemble system]
PostProcBrokenMeshInMoab< FaceElementForcesAndSourcesCore > PostProcFaceEle