v0.14.0
poisson_2d_nonhomogeneous.cpp
Go to the documentation of this file.
1 #include <stdlib.h>
4 
5 using namespace MoFEM;
7 
8 constexpr int SPACE_DIM = 2;
9 using PostProcFaceEle =
11 static char help[] = "...\n\n";
13 public:
15 
16  // Declaration of the main function to run analysis
17  MoFEMErrorCode runProgram();
18 
19 private:
20  // Declaration of other main functions called in runProgram()
21  MoFEMErrorCode readMesh();
22  MoFEMErrorCode setupProblem();
23  MoFEMErrorCode boundaryCondition();
24  MoFEMErrorCode assembleSystem();
25  MoFEMErrorCode setIntegrationRules();
26  MoFEMErrorCode solveSystem();
27  MoFEMErrorCode outputResults();
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]
62 
66 
68 }
69 //! [Read mesh]
70 
71 //! [Setup problem]
74 
79  int oRder = 3;
80  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &oRder, PETSC_NULL);
82 
84 
86 }
87 //! [Setup problem]
88 
89 //! [Boundary condition]
92 
93  auto bc_mng = mField.getInterface<BcManager>();
94  CHKERR bc_mng->pushMarkDOFsOnEntities(simpleInterface->getProblemName(),
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(
118  new OpSetBc(domainField, true, boundaryMarker));
119  pipeline_mng->getOpDomainLhsPipeline().push_back(
121  pipeline_mng->getOpDomainLhsPipeline().push_back(
122  new OpUnSetBc(domainField));
123  }
124 
125  { // Push operators to the Pipeline that is responsible for calculating RHS of
126  // domain elements
127  pipeline_mng->getOpDomainRhsPipeline().push_back(
128  new OpSetBc(domainField, true, boundaryMarker));
129  pipeline_mng->getOpDomainRhsPipeline().push_back(
131  pipeline_mng->getOpDomainRhsPipeline().push_back(
132  new OpUnSetBc(domainField));
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(
142  new OpUnSetBc(domainField));
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(
152  new OpUnSetBc(domainField));
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 
262 
263  readMesh();
264  setupProblem();
266  assembleSystem();
268  solveSystem();
269  outputResults();
270 
272 }
273 
274 int main(int argc, char *argv[]) {
275 
276  // Initialisation of MoFEM/PETSc and MOAB data structures
277  const char param_file[] = "param_file.petsc";
278  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
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  }
298  CATCH_ERRORS;
299 
300  // Finish work: cleaning memory, getting statistics, etc.
302 
303  return 0;
304 }
Poisson2DNonhomogeneousOperators::OpBoundaryRhs
[OpBoundaryLhs]
Definition: poisson_2d_nonhomogeneous.hpp:243
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
Poisson2DNonhomogeneousOperators::OpBoundaryLhs
[OpDomainRhs]
Definition: poisson_2d_nonhomogeneous.hpp:168
SPACE_DIM
constexpr int SPACE_DIM
Definition: poisson_2d_nonhomogeneous.cpp:8
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
Poisson2DNonhomogeneous::oRder
int oRder
Definition: poisson_2d_nonhomogeneous.cpp:48
main
int main(int argc, char *argv[])
Definition: poisson_2d_nonhomogeneous.cpp:274
Poisson2DNonhomogeneousOperators
Definition: poisson_2d_nonhomogeneous.hpp:15
Poisson2DNonhomogeneousOperators::OpDomainLhs
[OpDomainLhs]
Definition: poisson_2d_nonhomogeneous.hpp:22
Poisson2DNonhomogeneous
Definition: poisson_2d_nonhomogeneous.cpp:12
Poisson2DNonhomogeneousOperators::OpDomainRhs
[OpDomainLhs]
Definition: poisson_2d_nonhomogeneous.hpp:103
Poisson2DNonhomogeneous::simpleInterface
Simple * simpleInterface
Definition: poisson_2d_nonhomogeneous.cpp:44
Poisson2DNonhomogeneous::assembleSystem
MoFEMErrorCode assembleSystem()
[Boundary condition]
Definition: poisson_2d_nonhomogeneous.cpp:108
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
Poisson2DNonhomogeneous::mField
MoFEM::Interface & mField
Definition: poisson_2d_nonhomogeneous.cpp:43
MoFEM::Simple::loadFile
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
poisson_2d_nonhomogeneous.hpp
MoFEM::DMoFEMMeshToLocalVector
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:527
BasicFiniteElements.hpp
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Poisson2DNonhomogeneous::outputResults
MoFEMErrorCode outputResults()
[Solve system]
Definition: poisson_2d_nonhomogeneous.cpp:231
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
ROW
@ ROW
Definition: definitions.h:123
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:667
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1018
Poisson2DNonhomogeneous::Poisson2DNonhomogeneous
Poisson2DNonhomogeneous(MoFEM::Interface &m_field)
Definition: poisson_2d_nonhomogeneous.cpp:57
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
Poisson2DNonhomogeneous::readMesh
MoFEMErrorCode readMesh()
[Read mesh]
Definition: poisson_2d_nonhomogeneous.cpp:60
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
Poisson2DNonhomogeneous::domainField
std::string domainField
Definition: poisson_2d_nonhomogeneous.cpp:47
Poisson2DNonhomogeneous::boundaryFunction
static double boundaryFunction(const double x, const double y, const double z)
Definition: poisson_2d_nonhomogeneous.cpp:36
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
Poisson2DNonhomogeneous::sourceTermFunction
static double sourceTermFunction(const double x, const double y, const double z)
Definition: poisson_2d_nonhomogeneous.cpp:30
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::Simple::addDomainField
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
MoFEM::PipelineManager::createKSP
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
Definition: PipelineManager.cpp:49
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
Poisson2DNonhomogeneous::runProgram
MoFEMErrorCode runProgram()
Definition: poisson_2d_nonhomogeneous.cpp:260
Poisson2DNonhomogeneous::solveSystem
MoFEMErrorCode solveSystem()
[Set integration rules]
Definition: poisson_2d_nonhomogeneous.cpp:180
Poisson2DNonhomogeneous::setIntegrationRules
MoFEMErrorCode setIntegrationRules()
[Assemble system]
Definition: poisson_2d_nonhomogeneous.cpp:160
Poisson2DNonhomogeneous::boundaryEntitiesForFieldsplit
Range boundaryEntitiesForFieldsplit
Definition: poisson_2d_nonhomogeneous.cpp:54
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:473
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:503
MoFEM::OpUnSetBc
Definition: FormsIntegrators.hpp:49
Range
MoFEM::PipelineManager::getDomainLhsFE
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Definition: PipelineManager.hpp:401
MoFEM::CoreTmp< 0 >::Initialize
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
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:217
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1094
Poisson2DNonhomogeneous::boundaryCondition
MoFEMErrorCode boundaryCondition()
[Setup problem]
Definition: poisson_2d_nonhomogeneous.cpp:90
Poisson2DNonhomogeneous::setupProblem
MoFEMErrorCode setupProblem()
[Read mesh]
Definition: poisson_2d_nonhomogeneous.cpp:72
MoFEM::OpSetBc
Set indices on entities on finite element.
Definition: FormsIntegrators.hpp:38
MoFEM::DMMoFEMGetProblemPtr
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMoFEM.cpp:430
help
static char help[]
Definition: poisson_2d_nonhomogeneous.cpp:11
MoFEM::Simple::addBoundaryField
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
Poisson2DNonhomogeneous::boundaryMarker
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Definition: poisson_2d_nonhomogeneous.cpp:51
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:54
MoFEM::Simple::setUp
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:611
MoFEM::Problem::getName
auto getName() const
Definition: ProblemsMultiIndices.hpp:372
MoFEM::Simple::getProblemName
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:362
convert.int
int
Definition: convert.py:64
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
F
@ F
Definition: free_surface.cpp:394
MoFEM::PostProcBrokenMeshInMoab< FaceElementForcesAndSourcesCore >
Definition: PostProcBrokenMeshInMoabBase.hpp:677
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698