v0.14.0
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
15constexpr int SPACE_DIM =
16 EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
17
19
20using namespace MoFEM;
21using namespace Poisson2DHomogeneousOperators;
22
24
25static char help[] = "...\n\n";
26
28public:
30
31 // Declaration of the main function to run analysis
33
34private:
35 // Declaration of other main functions called in runProgram()
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 =
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
252
254}
255//! [Run program]
256
257//! [Main]
258int main(int argc, char *argv[]) {
259
260 // Initialisation of MoFEM/PETSc and MOAB data structures
261 const char param_file[] = "param_file.petsc";
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 }
283
284 // Finish work: cleaning memory, getting statistics, etc.
286
287 return 0;
288}
289//! [Main]
static Index< 'p', 3 > p
std::string param_file
int main()
Definition: adol-c_atom.cpp:46
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
static char help[]
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.
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
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:667
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:473
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:611
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:362
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]