v0.13.1
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
13static const int nb_ref_levels =
14 1; ///< if larger than zero set n-levels of random mesh refinements with
15 ///< hanging nodes
16
17constexpr auto field_name = "U";
18
22
23using namespace MoFEM;
24using namespace Poisson2DHomogeneousOperators;
25
28
29static char help[] = "...\n\n";
30
32public:
34
35 // Declaration of the main function to run analysis
37
38private:
39 // Declaration of other main functions called in runProgram()
47
48 // MoFEM interfaces
51
52 // Field name and approximation order
53 int oRder;
54};
55
57 : mField(m_field) {}
58
59//! [Read mesh]
62
66
68}
69//! [Read mesh]
70
71//! [Setup problem]
74
77
78 int oRder = 3;
79 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &oRder, PETSC_NULL);
81
82 // Refine random elements and create hanging nodes. This is only need if one
83 // like to refine mesh.
84 if (nb_ref_levels)
86
88
89 // Remove hanging nodes
92
94}
95//! [Setup problem]
96
97//! [Boundary condition]
100
101 // Get boundary edges marked in block named "BOUNDARY_CONDITION"
102 Range boundary_entities;
104 std::string entity_name = it->getName();
105 if (entity_name.compare(0, 18, "BOUNDARY_CONDITION") == 0) {
106 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 1,
107 boundary_entities, true);
108 }
109 }
110 // Add vertices to boundary entities
111 Range boundary_vertices;
112 CHKERR mField.get_moab().get_connectivity(boundary_entities,
113 boundary_vertices, true);
114 boundary_entities.merge(boundary_vertices);
115
116 // Remove DOFs as homogeneous boundary condition is used
117 CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
118 simpleInterface->getProblemName(), field_name, boundary_entities);
119
121}
122//! [Boundary condition]
123
124//! [Assemble system]
127
128 auto pipeline_mng = mField.getInterface<PipelineManager>();
129
130 auto det_ptr = boost::make_shared<VectorDouble>();
131 auto jac_ptr = boost::make_shared<MatrixDouble>();
132 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
133
134 { // Push operators to the Pipeline that is responsible for calculating LHS
135
136 pipeline_mng->getOpDomainLhsPipeline().push_back(
137 new OpCalculateHOJac<2>(jac_ptr));
138 pipeline_mng->getOpDomainLhsPipeline().push_back(
139 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
140 pipeline_mng->getOpDomainLhsPipeline().push_back(
141 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
142 pipeline_mng->getOpDomainLhsPipeline().push_back(
144
145 if (nb_ref_levels) { // This part is advanced. Can be skipped for not
146 // refined meshes with
147 // hanging nodes.
148 // Force integration on last refinement level, and add to top elements
149 // DOFs and based from underlying elements.
150 pipeline_mng->getDomainLhsFE()->exeTestHook = test_bit_child;
151 set_parent_dofs(mField, pipeline_mng->getDomainLhsFE(),
152 OpFaceEle::OPSPACE, QUIET, Sev::noisy);
153 set_parent_dofs(mField, pipeline_mng->getDomainLhsFE(), OpFaceEle::OPROW,
154 QUIET, Sev::noisy);
155 set_parent_dofs(mField, pipeline_mng->getDomainLhsFE(), OpFaceEle::OPCOL,
156 QUIET, Sev::noisy);
157 }
158
159 pipeline_mng->getOpDomainLhsPipeline().push_back(
161 }
162
163 { // Push operators to the Pipeline that is responsible for calculating RHS
164
165 if (nb_ref_levels) { // This part is advanced. Can be skipped for not
166 // refined meshes with
167 // hanging nodes.
168 // Force integration on last refinement level, and add to top elements
169 // DOFs and based from underlying elements.
170 pipeline_mng->getDomainRhsFE()->exeTestHook = test_bit_child;
171 set_parent_dofs(mField, pipeline_mng->getDomainRhsFE(),
172 OpFaceEle::OPSPACE, QUIET, Sev::noisy);
173 set_parent_dofs(mField, pipeline_mng->getDomainRhsFE(), OpFaceEle::OPROW,
174 QUIET, Sev::noisy);
175 }
176
177 pipeline_mng->getOpDomainRhsPipeline().push_back(
179 }
180
182}
183//! [Assemble system]
184
185//! [Set integration rules]
188
189 auto rule_lhs = [](int, int, int p) -> int { return 2 * (p - 1); };
190 auto rule_rhs = [](int, int, int p) -> int { return p; };
191
192 auto pipeline_mng = mField.getInterface<PipelineManager>();
193 CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule_lhs);
194 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule_rhs);
195
197}
198//! [Set integration rules]
199
200//! [Solve system]
203
204 auto pipeline_mng = mField.getInterface<PipelineManager>();
205
206 auto ksp_solver = pipeline_mng->createKSP();
207 CHKERR KSPSetFromOptions(ksp_solver);
208 CHKERR KSPSetUp(ksp_solver);
209
210 // Create RHS and solution vectors
211 auto dm = simpleInterface->getDM();
212 auto F = smartCreateDMVector(dm);
213 auto D = smartVectorDuplicate(F);
214
215 // Solve the system
216 CHKERR KSPSolve(ksp_solver, F, D);
217
218 // Scatter result data on the mesh
219 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
220 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
221 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
222
224}
225//! [Solve system]
226
227//! [Output results]
230
231 auto pipeline_mng = mField.getInterface<PipelineManager>();
232 pipeline_mng->getDomainLhsFE().reset();
233
234 auto post_proc_fe = boost::make_shared<PostProcFaceEle>(mField);
235
236 auto det_ptr = boost::make_shared<VectorDouble>();
237 auto jac_ptr = boost::make_shared<MatrixDouble>();
238 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
239
240 constexpr auto SPACE_DIM = 2; // dimension of problem
241
242 post_proc_fe->getOpPtrVector().push_back(
243 new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
244 post_proc_fe->getOpPtrVector().push_back(
245 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
246 post_proc_fe->getOpPtrVector().push_back(
247 new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
248
249 if (nb_ref_levels) { // This part is advanced. Can be skipped for not refined
250 // meshes with
251 // hanging nodes.
252 post_proc_fe->exeTestHook = test_bit_child;
254 Sev::noisy);
255 set_parent_dofs(mField, post_proc_fe, OpFaceEle::OPROW, QUIET, Sev::noisy);
256 }
257
258 auto u_ptr = boost::make_shared<VectorDouble>();
259 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
260 post_proc_fe->getOpPtrVector().push_back(
262
263 post_proc_fe->getOpPtrVector().push_back(
265
267
268 post_proc_fe->getOpPtrVector().push_back(
269
270 new OpPPMap(post_proc_fe->getPostProcMesh(),
271 post_proc_fe->getMapGaussPts(),
272
273 OpPPMap::DataMapVec{{"U", u_ptr}},
274
275 OpPPMap::DataMapMat{{"GRAD_U", grad_u_ptr}},
276
278
280
281 )
282
283 );
284
285 pipeline_mng->getDomainRhsFE() = post_proc_fe;
286 CHKERR pipeline_mng->loopFiniteElements();
287 CHKERR post_proc_fe->writeFile("out_result.h5m");
288
290}
291//! [Output results]
292
293//! [Run program]
296
304
306}
307//! [Run program]
308
309//! [Main]
310int main(int argc, char *argv[]) {
311
312 // Initialisation of MoFEM/PETSc and MOAB data structures
313 const char param_file[] = "param_file.petsc";
315
316 // Error handling
317 try {
318 // Register MoFEM discrete manager in PETSc
319 DMType dm_name = "DMMOFEM";
320 CHKERR DMRegister_MoFEM(dm_name);
321
322 // Create MOAB instance
323 moab::Core mb_instance; // mesh database
324 moab::Interface &moab = mb_instance; // mesh database interface
325
326 // Create MoFEM instance
327 MoFEM::Core core(moab); // finite element database
328 MoFEM::Interface &m_field = core; // finite element interface
329
330 // Run the main analysis
331 Poisson2DHomogeneous poisson_problem(m_field);
332 CHKERR poisson_problem.runProgram();
333 }
335
336 // Finish work: cleaning memory, getting statistics, etc.
338
339 return 0;
340}
341//! [Main]
static Index< 'p', 3 > p
std::string param_file
int main()
Definition: adol-c_atom.cpp:46
constexpr int SPACE_DIM
@ QUIET
Definition: definitions.h:208
#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
@ BLOCKSET
Definition: definitions.h:148
#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: DMMMoFEM.cpp:470
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
#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.
auto test_bit_child
lambda function used to select elements on which finite element pipelines are executed.
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.
auto set_parent_dofs
set levels of projection operators, which project field data from parent entities,...
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
static char help[]
constexpr auto field_name
static const int nb_ref_levels
virtual moab::Interface & get_moab()=0
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.
@ OPCOL
operator doWork function is executed on FE columns
@ OPROW
operator doWork function is executed on FE rows
@ OPSPACE
operator do Work is execute on space data
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.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
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:373
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:289
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:780
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:303
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:588
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:723
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:334
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]