v0.13.2
Loading...
Searching...
No Matches
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
13constexpr auto field_name = "U";
14
17
18using namespace MoFEM;
19using namespace Poisson2DHomogeneousOperators;
20
23
24static char help[] = "...\n\n";
25
27public:
29
30 // Declaration of the main function to run analysis
32
33private:
34 // Declaration of other main functions called in runProgram()
42
43 // MoFEM interfaces
46
47 // Field name and approximation order
48 int oRder;
49};
50
52 : mField(m_field) {}
53
54//! [Read mesh]
57
61
63}
64//! [Read mesh]
65
66//! [Setup problem]
69
72
73 int oRder = 3;
74 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &oRder, PETSC_NULL);
76
78
80}
81//! [Setup problem]
82
83//! [Boundary condition]
86
87 auto bc_mng = mField.getInterface<BcManager>();
88
89 // Remove BCs from blockset name "BOUNDARY_CONDITION";
90 CHKERR bc_mng->removeBlockDOFsOnEntities<BcScalarMeshsetType<BLOCKSET>>(
91 simpleInterface->getProblemName(), "BOUNDARY_CONDITION",
92 std::string(field_name), true);
93
95}
96//! [Boundary condition]
97
98//! [Assemble system]
101
102 auto pipeline_mng = mField.getInterface<PipelineManager>();
103
104 constexpr int space_dim = 2;
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 auto calculate_residual_from_set_values_on_bc = [&](auto &pipeline) {
124 using DomainEle =
129
130 auto grad_u_vals_ptr = boost::make_shared<MatrixDouble>();
131 pipeline_mng->getOpDomainRhsPipeline().push_back(
133 grad_u_vals_ptr));
134 pipeline_mng->getOpDomainRhsPipeline().push_back(
135 new OpInternal(field_name, grad_u_vals_ptr,
136 [](double, double, double) constexpr { return -1; }));
137 };
138
140 pipeline_mng->getOpDomainRhsPipeline(), {H1});
141 set_values_to_bc_dofs(pipeline_mng->getDomainRhsFE());
142 calculate_residual_from_set_values_on_bc(
143 pipeline_mng->getOpDomainRhsPipeline());
144 pipeline_mng->getOpDomainRhsPipeline().push_back(
146 }
147
149}
150//! [Assemble system]
151
152//! [Set integration rules]
155
156 auto rule_lhs = [](int, int, int p) -> int { return 2 * (p - 1); };
157 auto rule_rhs = [](int, int, int p) -> int { return p; };
158
159 auto pipeline_mng = mField.getInterface<PipelineManager>();
160 CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule_lhs);
161 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule_rhs);
162
164}
165//! [Set integration rules]
166
167//! [Solve system]
170
171 auto pipeline_mng = mField.getInterface<PipelineManager>();
172
173 auto ksp_solver = pipeline_mng->createKSP();
174 CHKERR KSPSetFromOptions(ksp_solver);
175 CHKERR KSPSetUp(ksp_solver);
176
177 // Create RHS and solution vectors
178 auto dm = simpleInterface->getDM();
179 auto F = smartCreateDMVector(dm);
180 auto D = smartVectorDuplicate(F);
181
182 // Solve the system
183 CHKERR KSPSolve(ksp_solver, F, D);
184
185 // Scatter result data on the mesh
186 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
187 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
188 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
189
191}
192//! [Solve system]
193
194//! [Output results]
197
198 auto pipeline_mng = mField.getInterface<PipelineManager>();
199 pipeline_mng->getDomainLhsFE().reset();
200
201 auto post_proc_fe = boost::make_shared<PostProcFaceEle>(mField);
202
203 auto det_ptr = boost::make_shared<VectorDouble>();
204 auto jac_ptr = boost::make_shared<MatrixDouble>();
205 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
206
207 constexpr auto SPACE_DIM = 2; // dimension of problem
208
209 post_proc_fe->getOpPtrVector().push_back(
210 new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
211 post_proc_fe->getOpPtrVector().push_back(
212 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
213 post_proc_fe->getOpPtrVector().push_back(
214 new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
215
216
217 auto u_ptr = boost::make_shared<VectorDouble>();
218 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
219 post_proc_fe->getOpPtrVector().push_back(
221
222 post_proc_fe->getOpPtrVector().push_back(
224
226
227 post_proc_fe->getOpPtrVector().push_back(
228
229 new OpPPMap(post_proc_fe->getPostProcMesh(),
230 post_proc_fe->getMapGaussPts(),
231
232 OpPPMap::DataMapVec{{"U", u_ptr}},
233
234 OpPPMap::DataMapMat{{"GRAD_U", grad_u_ptr}},
235
237
239
240 )
241
242 );
243
244 pipeline_mng->getDomainRhsFE() = post_proc_fe;
245 CHKERR pipeline_mng->loopFiniteElements();
246 CHKERR post_proc_fe->writeFile("out_result.h5m");
247
249}
250//! [Output results]
251
252//! [Run program]
255
263
265}
266//! [Run program]
267
268//! [Main]
269int main(int argc, char *argv[]) {
270
271 // Initialisation of MoFEM/PETSc and MOAB data structures
272 const char param_file[] = "param_file.petsc";
274
275 // Error handling
276 try {
277 // Register MoFEM discrete manager in PETSc
278 DMType dm_name = "DMMOFEM";
279 CHKERR DMRegister_MoFEM(dm_name);
280
281 // Create MOAB instance
282 moab::Core mb_instance; // mesh database
283 moab::Interface &moab = mb_instance; // mesh database interface
284
285 // Create MoFEM instance
286 MoFEM::Core core(moab); // finite element database
287 MoFEM::Interface &m_field = core; // finite element interface
288
289 // Run the main analysis
290 Poisson2DHomogeneous poisson_problem(m_field);
291 CHKERR poisson_problem.runProgram();
292 }
294
295 // Finish work: cleaning memory, getting statistics, etc.
297
298 return 0;
299}
300//! [Main]
static Index< 'p', 3 > p
std::string param_file
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
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:511
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:995
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 > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
static char help[]
constexpr auto field_name
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition: BcManager.hpp:24
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:39
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.
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
Set inverse jacobian to base functions.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
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
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:671
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:194
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:476
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:614
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:348
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]