v0.14.0
poisson_2d_homogeneous.cpp
Go to the documentation of this file.
1 /**
2  * \file poisson_2d_homogeneous.cpp
3  * \example poisson_2d_homogeneous.cpp
4  *
5  * Solution of poisson equation. Direct implementation of User Data Operators
6  * for teaching proposes.
7  *
8  * \note In practical application we suggest use form integrators to generalise
9  * and simplify code. However, here we like to expose user to ways how to
10  * implement data operator from scratch.
11  */
12 
13 constexpr auto field_name = "U";
14 
15 constexpr int SPACE_DIM =
16  EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
17 
19 
20 using namespace MoFEM;
21 using namespace Poisson2DHomogeneousOperators;
22 
24 
25 static char help[] = "...\n\n";
26 
28 public:
30 
31  // Declaration of the main function to run analysis
32  MoFEMErrorCode runProgram();
33 
34 private:
35  // Declaration of other main functions called in runProgram()
36  MoFEMErrorCode readMesh();
37  MoFEMErrorCode setupProblem();
38  MoFEMErrorCode boundaryCondition();
39  MoFEMErrorCode assembleSystem();
40  MoFEMErrorCode setIntegrationRules();
41  MoFEMErrorCode solveSystem();
42  MoFEMErrorCode outputResults();
43 
44  // MoFEM interfaces
47 
48  // Field name and approximation order
49  int oRder;
50 };
51 
53  : mField(m_field) {}
54 
55 //! [Read mesh]
58 
62 
64 }
65 //! [Read mesh]
66 
67 //! [Setup problem]
70 
73 
74  int oRder = 3;
75  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &oRder, PETSC_NULL);
77 
79 
81 }
82 //! [Setup problem]
83 
84 //! [Boundary condition]
87 
88  auto bc_mng = mField.getInterface<BcManager>();
89 
90  // Remove BCs from blockset name "BOUNDARY_CONDITION" or SETU, note that you
91  // can use regular expression to put list of blocksets;
92  CHKERR bc_mng->removeBlockDOFsOnEntities<BcScalarMeshsetType<BLOCKSET>>(
93  simpleInterface->getProblemName(), "(BOUNDARY_CONDITION|SETU)",
94  std::string(field_name), true);
95 
97 }
98 //! [Boundary condition]
99 
100 //! [Assemble system]
103 
104  auto pipeline_mng = mField.getInterface<PipelineManager>();
105 
106  { // Push operators to the Pipeline that is responsible for calculating LHS
108  pipeline_mng->getOpDomainLhsPipeline(), {H1});
109  pipeline_mng->getOpDomainLhsPipeline().push_back(
111  }
112 
113  { // Push operators to the Pipeline that is responsible for calculating RHS
114 
115  auto set_values_to_bc_dofs = [&](auto &fe) {
116  auto get_bc_hook = [&]() {
118  return hook;
119  };
120  fe->preProcessHook = get_bc_hook();
121  };
122 
123  // you can skip that if boundary condition is prescribing zero
124  auto calculate_residual_from_set_values_on_bc = [&](auto &pipeline) {
125  using DomainEle =
128  using OpInternal = FormsIntegrators<DomainEleOp>::Assembly<
130 
131  auto grad_u_vals_ptr = boost::make_shared<MatrixDouble>();
132  pipeline_mng->getOpDomainRhsPipeline().push_back(
134  grad_u_vals_ptr));
135  pipeline_mng->getOpDomainRhsPipeline().push_back(
136  new OpInternal(field_name, grad_u_vals_ptr,
137  [](double, double, double) constexpr { return -1; }));
138  };
139 
141  pipeline_mng->getOpDomainRhsPipeline(), {H1});
142  set_values_to_bc_dofs(pipeline_mng->getDomainRhsFE());
143  calculate_residual_from_set_values_on_bc(
144  pipeline_mng->getOpDomainRhsPipeline());
145  pipeline_mng->getOpDomainRhsPipeline().push_back(
147  }
148 
150 }
151 //! [Assemble system]
152 
153 //! [Set integration rules]
156 
157  auto rule_lhs = [](int, int, int p) -> int { return 2 * (p - 1); };
158  auto rule_rhs = [](int, int, int p) -> int { return p; };
159 
160  auto pipeline_mng = mField.getInterface<PipelineManager>();
161  CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule_lhs);
162  CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule_rhs);
163 
165 }
166 //! [Set integration rules]
167 
168 //! [Solve system]
171 
172  auto pipeline_mng = mField.getInterface<PipelineManager>();
173 
174  auto ksp_solver = pipeline_mng->createKSP();
175  CHKERR KSPSetFromOptions(ksp_solver);
176  CHKERR KSPSetUp(ksp_solver);
177 
178  // Create RHS and solution vectors
179  auto dm = simpleInterface->getDM();
180  auto F = createDMVector(dm);
181  auto D = vectorDuplicate(F);
182 
183  // Solve the system
184  CHKERR KSPSolve(ksp_solver, F, D);
185 
186  // Scatter result data on the mesh
187  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
188  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
189  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
190 
192 }
193 //! [Solve system]
194 
195 //! [Output results]
198 
199  auto pipeline_mng = mField.getInterface<PipelineManager>();
200  pipeline_mng->getDomainLhsFE().reset();
201 
202  auto post_proc_fe = boost::make_shared<PostProcFaceEle>(mField);
204  post_proc_fe->getOpPtrVector(), {H1});
205 
206  auto u_ptr = boost::make_shared<VectorDouble>();
207  auto grad_u_ptr = boost::make_shared<MatrixDouble>();
208  post_proc_fe->getOpPtrVector().push_back(
210 
211  post_proc_fe->getOpPtrVector().push_back(
213 
215 
216  post_proc_fe->getOpPtrVector().push_back(
217 
218  new OpPPMap(post_proc_fe->getPostProcMesh(),
219  post_proc_fe->getMapGaussPts(),
220 
221  OpPPMap::DataMapVec{{"U", u_ptr}},
222 
223  OpPPMap::DataMapMat{{"GRAD_U", grad_u_ptr}},
224 
226 
228 
229  )
230 
231  );
232 
233  pipeline_mng->getDomainRhsFE() = post_proc_fe;
234  CHKERR pipeline_mng->loopFiniteElements();
235  CHKERR post_proc_fe->writeFile("out_result.h5m");
236 
238 }
239 //! [Output results]
240 
241 //! [Run program]
244 
245  CHKERR readMesh();
252 
254 }
255 //! [Run program]
256 
257 //! [Main]
258 int main(int argc, char *argv[]) {
259 
260  // Initialisation of MoFEM/PETSc and MOAB data structures
261  const char param_file[] = "param_file.petsc";
262  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
263 
264  // Error handling
265  try {
266  // Register MoFEM discrete manager in PETSc
267  DMType dm_name = "DMMOFEM";
268  CHKERR DMRegister_MoFEM(dm_name);
269 
270  // Create MOAB instance
271  moab::Core mb_instance; // mesh database
272  moab::Interface &moab = mb_instance; // mesh database interface
273 
274  // Create MoFEM instance
275  MoFEM::Core core(moab); // finite element database
276  MoFEM::Interface &m_field = core; // finite element interface
277 
278  // Run the main analysis
279  Poisson2DHomogeneous poisson_problem(m_field);
280  CHKERR poisson_problem.runProgram();
281  }
282  CATCH_ERRORS;
283 
284  // Finish work: cleaning memory, getting statistics, etc.
286 
287  return 0;
288 }
289 //! [Main]
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
Poisson2DHomogeneous::solveSystem
MoFEMErrorCode solveSystem()
[Set integration rules]
Definition: poisson_2d_homogeneous.cpp:169
EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
Poisson2DHomogeneousOperators::OpDomainRhsVectorF
Definition: poisson_2d_homogeneous.hpp:92
Poisson2DHomogeneous::oRder
int oRder
Definition: poisson_2d_homogeneous.cpp:49
MoFEM::PipelineManager::ElementsAndOpsByDim
Definition: PipelineManager.hpp:38
Poisson2DHomogeneous::setupProblem
MoFEMErrorCode setupProblem()
[Read mesh]
Definition: poisson_2d_homogeneous.cpp:68
MoFEM::OpPostProcMapInMoab::DataMapVec
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
Definition: PostProcBrokenMeshInMoabBase.hpp:700
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Simple::loadFile
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
Poisson2DHomogeneous::assembleSystem
MoFEMErrorCode assembleSystem()
[Boundary condition]
Definition: poisson_2d_homogeneous.cpp:101
Poisson2DHomogeneous::outputResults
MoFEMErrorCode outputResults()
[Solve system]
Definition: poisson_2d_homogeneous.cpp:196
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
SPACE_DIM
constexpr int SPACE_DIM
Definition: poisson_2d_homogeneous.cpp:15
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Poisson2DHomogeneous::Poisson2DHomogeneous
Poisson2DHomogeneous(MoFEM::Interface &m_field)
Definition: poisson_2d_homogeneous.cpp:52
Poisson2DHomogeneous::simpleInterface
Simple * simpleInterface
Definition: poisson_2d_homogeneous.cpp:46
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1294
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:667
Poisson2DHomogeneousOperators::OpDomainLhsMatrixK
Definition: poisson_2d_homogeneous.hpp:38
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:667
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1067
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:310
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::BcScalarMeshsetType
Definition: BcManager.hpp:15
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::Simple::addDomainField
MoFEMErrorCode addDomainField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:264
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
Poisson2DHomogeneous
Definition: poisson_2d_homogeneous.cpp:27
Poisson2DHomogeneous::mField
MoFEM::Interface & mField
Definition: poisson_2d_homogeneous.cpp:45
Poisson2DHomogeneous::setIntegrationRules
MoFEMErrorCode setIntegrationRules()
[Assemble system]
Definition: poisson_2d_homogeneous.cpp:154
MoFEM::OpPostProcMapInMoab::DataMapMat
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
Definition: PostProcBrokenMeshInMoabBase.hpp:701
main
int main(int argc, char *argv[])
[Run program]
Definition: poisson_2d_homogeneous.cpp:258
OpGradTimesTensor
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
Definition: operators_tests.cpp:48
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
help
static char help[]
Definition: poisson_2d_homogeneous.cpp:25
Poisson2DHomogeneous::boundaryCondition
MoFEMErrorCode boundaryCondition()
[Setup problem]
Definition: poisson_2d_homogeneous.cpp:85
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:473
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:503
MoFEM::EssentialPreProc< TemperatureCubitBcData >
Specialization for TemperatureCubitBcData.
Definition: EssentialTemperatureCubitBcData.hpp:91
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::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
Poisson2DHomogeneous::runProgram
MoFEMErrorCode runProgram()
[Output results]
Definition: poisson_2d_homogeneous.cpp:242
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1102
Poisson2DHomogeneousOperators
Definition: poisson_2d_homogeneous.hpp:30
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
poisson_2d_homogeneous.hpp
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
Poisson2DHomogeneous::readMesh
MoFEMErrorCode readMesh()
[Read mesh]
Definition: poisson_2d_homogeneous.cpp:56
MoFEM::Simple::setUp
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:611
MoFEM::Simple::getProblemName
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:362
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:416
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
F
@ F
Definition: free_surface.cpp:394
MoFEM::PostProcBrokenMeshInMoab< FaceElementForcesAndSourcesCore >
Definition: PostProcBrokenMeshInMoabBase.hpp:677
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698