v0.13.1
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
19#include <BasicFiniteElements.hpp>
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};
56
58 : mField(m_field) {}
59
60//! [Read mesh]
63
67
69}
70//! [Read mesh]
71
72//! [Setup problem]
75
78
79 int oRder = 3;
80 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &oRder, PETSC_NULL);
82
83 // Refine random elements and create hanging nodes. This is only need if one
84 // like to refine mesh.
85 if (nb_ref_levels)
87
89
90 // Remove hanging nodes
93
95}
96//! [Setup problem]
97
98//! [Boundary condition]
101
102 // Get boundary edges marked in block named "BOUNDARY_CONDITION"
103 Range boundary_entities;
105 std::string entity_name = it->getName();
106 if (entity_name.compare(0, 18, "BOUNDARY_CONDITION") == 0) {
107 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 1,
108 boundary_entities, true);
109 }
110 }
111 // Add vertices to boundary entities
112 Range boundary_vertices;
113 CHKERR mField.get_moab().get_connectivity(boundary_entities,
114 boundary_vertices, true);
115 boundary_entities.merge(boundary_vertices);
116
117 // Remove DOFs as homogeneous boundary condition is used
118 CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
119 simpleInterface->getProblemName(), field_name, boundary_entities);
120
122}
123//! [Boundary condition]
124
125//! [Assemble system]
128
129 auto pipeline_mng = mField.getInterface<PipelineManager>();
130
131 auto det_ptr = boost::make_shared<VectorDouble>();
132 auto jac_ptr = boost::make_shared<MatrixDouble>();
133 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
134
135 { // Push operators to the Pipeline that is responsible for calculating LHS
136
137 pipeline_mng->getOpDomainLhsPipeline().push_back(
138 new OpCalculateHOJac<2>(jac_ptr));
139 pipeline_mng->getOpDomainLhsPipeline().push_back(
140 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
141 pipeline_mng->getOpDomainLhsPipeline().push_back(
142 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
143 pipeline_mng->getOpDomainLhsPipeline().push_back(
145
146 if (nb_ref_levels) { // This part is advanced. Can be skipped for not
147 // refined meshes with
148 // hanging nodes.
149 // Force integration on last refinement level, and add to top elements
150 // DOFs and based from underlying elements.
151 pipeline_mng->getDomainLhsFE()->exeTestHook = test_bit_child;
152 set_parent_dofs(mField, pipeline_mng->getDomainLhsFE(),
153 OpFaceEle::OPSPACE, QUIET, Sev::noisy);
154 set_parent_dofs(mField, pipeline_mng->getDomainLhsFE(), OpFaceEle::OPROW,
155 QUIET, Sev::noisy);
156 set_parent_dofs(mField, pipeline_mng->getDomainLhsFE(), OpFaceEle::OPCOL,
157 QUIET, Sev::noisy);
158 }
159
160 pipeline_mng->getOpDomainLhsPipeline().push_back(
162 }
163
164 { // Push operators to the Pipeline that is responsible for calculating RHS
165
166 if (nb_ref_levels) { // This part is advanced. Can be skipped for not
167 // refined meshes with
168 // hanging nodes.
169 // Force integration on last refinement level, and add to top elements
170 // DOFs and based from underlying elements.
171 pipeline_mng->getDomainRhsFE()->exeTestHook = test_bit_child;
172 set_parent_dofs(mField, pipeline_mng->getDomainRhsFE(),
173 OpFaceEle::OPSPACE, QUIET, Sev::noisy);
174 set_parent_dofs(mField, pipeline_mng->getDomainRhsFE(), OpFaceEle::OPROW,
175 QUIET, Sev::noisy);
176 }
177
178 pipeline_mng->getOpDomainRhsPipeline().push_back(
180 }
181
183}
184//! [Assemble system]
185
186//! [Set integration rules]
189
190 auto rule_lhs = [](int, int, int p) -> int { return 2 * (p - 1); };
191 auto rule_rhs = [](int, int, int p) -> int { return p; };
192
193 auto pipeline_mng = mField.getInterface<PipelineManager>();
194 CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule_lhs);
195 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule_rhs);
196
198}
199//! [Set integration rules]
200
201//! [Solve system]
204
205 auto pipeline_mng = mField.getInterface<PipelineManager>();
206
207 auto ksp_solver = pipeline_mng->createKSP();
208 CHKERR KSPSetFromOptions(ksp_solver);
209 CHKERR KSPSetUp(ksp_solver);
210
211 // Create RHS and solution vectors
212 auto dm = simpleInterface->getDM();
213 auto F = smartCreateDMVector(dm);
214 auto D = smartVectorDuplicate(F);
215
216 // Solve the system
217 CHKERR KSPSolve(ksp_solver, F, D);
218
219 // Scatter result data on the mesh
220 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
221 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
222 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
223
225}
226//! [Solve system]
227
228//! [Output results]
231
232 auto pipeline_mng = mField.getInterface<PipelineManager>();
233 pipeline_mng->getDomainLhsFE().reset();
234
235 auto post_proc_fe = boost::make_shared<PostProcFaceEle>(mField);
236
237 auto det_ptr = boost::make_shared<VectorDouble>();
238 auto jac_ptr = boost::make_shared<MatrixDouble>();
239 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
240
241 constexpr auto SPACE_DIM = 2; // dimension of problem
242
243 post_proc_fe->getOpPtrVector().push_back(
244 new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
245 post_proc_fe->getOpPtrVector().push_back(
246 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
247 post_proc_fe->getOpPtrVector().push_back(
248 new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
249
250 if (nb_ref_levels) { // This part is advanced. Can be skipped for not refined
251 // meshes with
252 // hanging nodes.
253 post_proc_fe->exeTestHook = test_bit_child;
255 Sev::noisy);
256 set_parent_dofs(mField, post_proc_fe, OpFaceEle::OPROW, QUIET, Sev::noisy);
257 }
258
259 auto u_ptr = boost::make_shared<VectorDouble>();
260 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
261 post_proc_fe->getOpPtrVector().push_back(
263
264 post_proc_fe->getOpPtrVector().push_back(
266
268
269 post_proc_fe->getOpPtrVector().push_back(
270
271 new OpPPMap(post_proc_fe->getPostProcMesh(),
272 post_proc_fe->getMapGaussPts(),
273
274 OpPPMap::DataMapVec{{"U", u_ptr}},
275
276 OpPPMap::DataMapMat{{"GRAD_U", grad_u_ptr}},
277
279
281
282 )
283
284 );
285
286 pipeline_mng->getDomainRhsFE() = post_proc_fe;
287 CHKERR pipeline_mng->loopFiniteElements();
288 CHKERR post_proc_fe->writeFile("out_result.h5m");
289
291}
292//! [Output results]
293
294//! [Run program]
297
305
307}
308//! [Run program]
309
310//! [Main]
311int main(int argc, char *argv[]) {
312
313 // Initialisation of MoFEM/PETSc and MOAB data structures
314 const char param_file[] = "param_file.petsc";
316
317 // Error handling
318 try {
319 // Register MoFEM discrete manager in PETSc
320 DMType dm_name = "DMMOFEM";
321 CHKERR DMRegister_MoFEM(dm_name);
322
323 // Create MOAB instance
324 moab::Core mb_instance; // mesh database
325 moab::Interface &moab = mb_instance; // mesh database interface
326
327 // Create MoFEM instance
328 MoFEM::Core core(moab); // finite element database
329 MoFEM::Interface &m_field = core; // finite element interface
330
331 // Run the main analysis
332 Poisson2DHomogeneous poisson_problem(m_field);
333 CHKERR poisson_problem.runProgram();
334 }
336
337 // Finish work: cleaning memory, getting statistics, etc.
339
340 return 0;
341}
342//! [Main]
static Index< 'p', 3 > p
std::string param_file
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.
auto set_parent_dofs(MoFEM::Interface &m_field, boost::shared_ptr< FEMethod > &fe_top, ForcesAndSourcesCore::UserDataOperator::OpType op, int verbosity, LogManager::SeverityLevel sev)
set levels of projection operators, which project field data from parent entities,...
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
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.
CoreTmp< 0 > Core
Definition: Core.hpp:1086
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
const double D
diffusivity
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
int main(int argc, char *argv[])
[Run program]
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.
Definition: PostProc.hpp:166
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
Definition: PostProc.hpp:168
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
Definition: PostProc.hpp:169
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()
Definition: PostProc.hpp:131
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:374
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:290
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:806
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:304
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:589
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:749
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]