v0.14.0
Loading...
Searching...
No Matches
poisson_2d_nonhomogeneous.cpp
Go to the documentation of this file.
1#include <stdlib.h>
4
5using namespace MoFEM;
7
8constexpr int SPACE_DIM = 2;
11static char help[] = "...\n\n";
13public:
15
16 // Declaration of the main function to run analysis
18
19private:
20 // Declaration of other main functions called in runProgram()
28
29 // Function to calculate the Source term
30 static double sourceTermFunction(const double x, const double y,
31 const double z) {
32 return 200 * sin(x * 10.) * cos(y * 10.);
33 // return 1;
34 }
35 // Function to calculate the Boundary term
36 static double boundaryFunction(const double x, const double y,
37 const double z) {
38 return sin(x * 10.) * cos(y * 10.);
39 // return 0;
40 }
41
42 // Main interfaces
45
46 // Field name and approximation order
47 std::string domainField;
48 int oRder;
49
50 // Object to mark boundary entities for the assembling of domain elements
51 boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
52
53 // Boundary entities marked for fieldsplit (block) solver - optional
55};
56
58 : domainField("U"), mField(m_field) {}
59//! [Read mesh]
69//! [Read mesh]
70
71//! [Setup problem]
87//! [Setup problem]
88
89//! [Boundary condition]
92
93 auto bc_mng = mField.getInterface<BcManager>();
95 "BOUNDARY_CONDITION", domainField, 0, 1,
96 true);
97 // merge markers from all blocksets "BOUNDARY_CONDITION"
98 boundaryMarker = bc_mng->getMergedBlocksMarker({"BOUNDARY_CONDITION"});
99 // get entities on blocksets "BOUNDARY_CONDITION"
101 bc_mng->getMergedBlocksRange({"BOUNDARY_CONDITION"});
102
104}
105//! [Boundary condition]
106
107//! [Assemble system]
110
111 auto pipeline_mng = mField.getInterface<PipelineManager>();
112
113 { // Push operators to the Pipeline that is responsible for calculating LHS of
114 // domain elements
116 pipeline_mng->getOpDomainLhsPipeline(), {H1});
117 pipeline_mng->getOpDomainLhsPipeline().push_back(
119 pipeline_mng->getOpDomainLhsPipeline().push_back(
121 pipeline_mng->getOpDomainLhsPipeline().push_back(
123 }
124
125 { // Push operators to the Pipeline that is responsible for calculating RHS of
126 // domain elements
127 pipeline_mng->getOpDomainRhsPipeline().push_back(
129 pipeline_mng->getOpDomainRhsPipeline().push_back(
131 pipeline_mng->getOpDomainRhsPipeline().push_back(
133 }
134
135 { // Push operators to the Pipeline that is responsible for calculating LHS of
136 // boundary elements
137 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
138 new OpSetBc(domainField, false, boundaryMarker));
139 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
141 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
143 }
144
145 { // Push operators to the Pipeline that is responsible for calculating RHS of
146 // boundary elements
147 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
148 new OpSetBc(domainField, false, boundaryMarker));
149 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
151 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
153 }
154
156}
157//! [Assemble system]
158
159//! [Set integration rules]
162
163 auto pipeline_mng = mField.getInterface<PipelineManager>();
164
165 auto domain_rule_lhs = [](int, int, int p) -> int { return 2 * (p - 1); };
166 auto domain_rule_rhs = [](int, int, int p) -> int { return 2 * (p - 1); };
167 CHKERR pipeline_mng->setDomainLhsIntegrationRule(domain_rule_lhs);
168 CHKERR pipeline_mng->setDomainRhsIntegrationRule(domain_rule_rhs);
169
170 auto boundary_rule_lhs = [](int, int, int p) -> int { return 2 * p; };
171 auto boundary_rule_rhs = [](int, int, int p) -> int { return 2 * p; };
172 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(boundary_rule_lhs);
173 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(boundary_rule_rhs);
174
176}
177//! [Set integration rules]
178
179//! [Solve system]
182
183 auto pipeline_mng = mField.getInterface<PipelineManager>();
184
185 auto ksp_solver = pipeline_mng->createKSP();
186 CHKERR KSPSetFromOptions(ksp_solver);
187
188 // Create RHS and solution vectors
189 auto dm = simpleInterface->getDM();
190 auto F = createDMVector(dm);
191 auto D = vectorDuplicate(F);
192
193 // Setup fieldsplit (block) solver - optional: yes/no
194 if (1) {
195 PC pc;
196 CHKERR KSPGetPC(ksp_solver, &pc);
197 PetscBool is_pcfs = PETSC_FALSE;
198 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
199
200 // Set up FIELDSPLIT, only when user set -pc_type fieldsplit
201 // Identify the index for boundary entities, remaining will be for domain
202 // Then split the fields for boundary and domain for solving
203 if (is_pcfs == PETSC_TRUE) {
204 IS is_domain, is_boundary;
205 const MoFEM::Problem *problem_ptr;
206 CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
207 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
208 problem_ptr->getName(), ROW, domainField, 0, 1, &is_boundary,
210 // CHKERR ISView(is_boundary, PETSC_VIEWER_STDOUT_SELF);
211 CHKERR PCFieldSplitSetIS(pc, NULL, is_boundary);
212 CHKERR ISDestroy(&is_boundary);
213 }
214 }
215
216 CHKERR KSPSetUp(ksp_solver);
217
218 // Solve the system
219 CHKERR KSPSolve(ksp_solver, F, D);
220
221 // Scatter result data on the mesh
222 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
223 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
224 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
225
227}
228//! [Solve system]
229
230//! [Output results]
233
234 auto pipeline_mng = mField.getInterface<PipelineManager>();
235 pipeline_mng->getDomainLhsFE().reset();
236 pipeline_mng->getBoundaryLhsFE().reset();
237 pipeline_mng->getDomainRhsFE().reset();
238 pipeline_mng->getBoundaryRhsFE().reset();
239
240 auto d_ptr = boost::make_shared<VectorDouble>();
241 auto l_ptr = boost::make_shared<VectorDouble>();
242
244
245 auto post_proc_fe = boost::make_shared<PostProcFaceEle>(mField);
246
247 post_proc_fe->getOpPtrVector().push_back(
249 post_proc_fe->getOpPtrVector().push_back(new OpPPMap(
250 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
251 {{domainField, d_ptr}}, {}, {}, {}));
252 pipeline_mng->getDomainRhsFE() = post_proc_fe;
253
254 CHKERR pipeline_mng->loopFiniteElements();
255 CHKERR post_proc_fe->writeFile("out_result.h5m");
256
258}
259
273
274int main(int argc, char *argv[]) {
275
276 // Initialisation of MoFEM/PETSc and MOAB data structures
277 const char param_file[] = "param_file.petsc";
279
280 // Error handling
281 try {
282 // Register MoFEM discrete manager in PETSc
283 DMType dm_name = "DMMOFEM";
284 CHKERR DMRegister_MoFEM(dm_name);
285
286 // Create MOAB instance
287 moab::Core mb_instance; // mesh database
288 moab::Interface &moab = mb_instance; // mesh database interface
289
290 // Create MoFEM instance
291 MoFEM::Core core(moab); // finite element database
292 MoFEM::Interface &m_field = core; // finite element interface
293
294 // Run the main analysis
295 Poisson2DNonhomogeneous poisson_problem(m_field);
296 CHKERR poisson_problem.runProgram();
297 }
299
300 // Finish work: cleaning memory, getting statistics, etc.
302
303 return 0;
304}
static Index< 'p', 3 > p
std::string param_file
int main()
@ ROW
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition definitions.h:64
@ 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 ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
@ 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 DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition DMMoFEM.cpp:412
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.
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
static char help[]
constexpr int SPACE_DIM
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition BcManager.hpp:25
MoFEMErrorCode pushMarkDOFsOnEntities(const std::string problem_name, const std::string block_name, const std::string field_name, int lo, int hi, bool get_low_dim_ents=true)
Mark block DOFs.
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.
Section manager is used to create indexes and sections.
Definition ISManager.hpp:23
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
Set indices on entities on finite element.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
keeps basic data about problem
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 addBoundaryField(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 boundary.
Definition Simple.cpp:282
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.
static double boundaryFunction(const double x, const double y, const double z)
Poisson2DNonhomogeneous(MoFEM::Interface &m_field)
MoFEMErrorCode solveSystem()
[Set integration rules]
MoFEMErrorCode setIntegrationRules()
[Assemble system]
static double sourceTermFunction(const double x, const double y, const double z)
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode boundaryCondition()
[Setup problem]
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
MoFEMErrorCode outputResults()
[Solve system]
MoFEMErrorCode readMesh()
[Read mesh]
MoFEMErrorCode assembleSystem()
[Boundary condition]