v0.13.1
Loading...
Searching...
No Matches
poisson_prob.cpp
Go to the documentation of this file.
1#include <stdlib.h>
3#include <PoissonOperators.hpp>
4
5using namespace MoFEM;
6using namespace PoissonOps;
7
8static char help[] = "...\n\n";
9
10struct Poisson {
11public:
12 Poisson(moab::Core &mb_instance, MoFEM::Core &core, int order);
13
15
16private:
24
25 boost::shared_ptr<VolEle> volPipelineLhs;
26 boost::shared_ptr<VolEle> volPipelineRhs;
27 boost::shared_ptr<VolEle> postFaceVolume;
28
30 moab::Interface &moab;
32
33 DM dm;
34 KSP ksp;
35
37
38 MPI_Comm mpiComm; // mpi parallel communicator
39 const int mpiRank; // number of processors involved
40
41 int order;
42
43 std::string domainField;
44};
45
46Poisson::Poisson(moab::Core &mb_instance, MoFEM::Core &core, const int order)
47 : domainField("U"), moab(mb_instance), mField(core),
48 mpiComm(mField.get_comm()), mpiRank(mField.get_comm_rank()),
49 order(order) {
50 volPipelineLhs = boost::shared_ptr<VolEle>(new VolEle(mField));
51 volPipelineRhs = boost::shared_ptr<VolEle>(new VolEle(mField));
52}
53
56
60
62}
63
66
71
73}
74
77 auto vol_ruleLhs = [](int, int, int p) -> int { return 2 * (p - 1); };
78 auto vol_ruleRhs = [](int, int, int p) -> int { return 2 * (p - 1); };
79 volPipelineLhs->getRuleHook = vol_ruleLhs;
80 volPipelineRhs->getRuleHook = vol_ruleRhs;
82}
83
86 auto det_ptr = boost::make_shared<VectorDouble>();
87 auto jac_ptr = boost::make_shared<MatrixDouble>();
88 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
89 volPipelineLhs->getOpPtrVector().push_back(
90 new OpCalculateHOJacForFace(jac_ptr));
91 volPipelineLhs->getOpPtrVector().push_back(
92 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
93 volPipelineLhs->getOpPtrVector().push_back(
95 volPipelineRhs->getOpPtrVector().push_back(new OpRhsU(domainField));
96
97 // dm = simpleInterface->getDM();
99 boost::shared_ptr<VolEle> null;
100
102
106}
107
110
111 Vec global_rhs, global_solution;
112 CHKERR DMCreateGlobalVector(dm, &global_rhs);
113 CHKERR VecDuplicate(global_rhs, &global_solution);
114
115
116 // ksp = createKSP(mField.get_comm());
117
118 CHKERR KSPCreate(PETSC_COMM_WORLD, &ksp);
119
120 // CHKERR KSPSetOperators(ksp, volPipelineLhs->getFEMethod()->ksp_B,
121 // volPipelineLhs->getFEMethod()->ksp_B);
122
123 CHKERR KSPSetFromOptions(ksp);
124
125 CHKERR KSPSetDM(ksp, dm);
126
127 CHKERR KSPSetUp(ksp);
128
129 CHKERR KSPSolve(ksp, global_rhs, global_solution);
130
131 CHKERR DMoFEMMeshToGlobalVector(dm, global_solution, INSERT_VALUES,
132 SCATTER_REVERSE);
133
134
135
136 CHKERR KSPDestroy(&ksp);
137 CHKERR VecDestroy(&global_solution);
138 CHKERR VecDestroy(&global_rhs);
139
142
143 CHKERR
144 boost::static_pointer_cast<PostProcFaceOnRefinedMesh>(postFaceVolume)
145 ->writeFile("out_vol.h5m");
146
147 CHKERR DMDestroy(&dm);
148
150}
151
154
156 boost::shared_ptr<VolEle>(new PostProcFaceOnRefinedMesh(mField));
157
158 CHKERR
159 boost::static_pointer_cast<PostProcFaceOnRefinedMesh>(postFaceVolume)
160 ->generateReferenceElementMesh();
161 CHKERR
162 boost::static_pointer_cast<PostProcFaceOnRefinedMesh>(postFaceVolume)
163 ->addFieldValuesPostProc(domainField);
164 // CHKERR
165 // boost::static_pointer_cast<PostProcFaceOnRefinedMesh>(postFaceVolume)
166 // ->addFieldValuesGradientPostProc(domainField);
167
169}
170
173 Range bdry_entities;
175 string name = it->getName();
176 if (name.compare(0, 14, "ESSENTIAL") == 0) {
177 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 1,
178 bdry_entities, true);
179 }
180 }
181
182 cerr << bdry_entities;
183 Range bdr_verts;
184 CHKERR mField.get_moab().get_connectivity(bdry_entities, bdr_verts, true);
185 bdry_entities.merge(bdr_verts);
186
187 CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
188 simpleInterface->getProblemName(), domainField, bdry_entities);
189
190
192}
193
196
199 CHKERR applyBc();
200
202
204
206
208}
209
210int main(int argc, char *argv[]) {
211 const char param_file[] = "param_file.petsc";
213 try {
214 moab::Core mb_instance;
215 moab::Interface &moab = mb_instance;
216 MoFEM::Core core(moab);
217 DMType dm_name = "DMMOFEM";
218 CHKERR DMRegister_MoFEM(dm_name);
219
220 int order = 1;
221 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
222
223 Poisson poisson_problem(mb_instance, core, order);
224
225 CHKERR poisson_problem.runAnalysis();
226 }
229 return 0;
230}
static Index< 'p', 3 > p
std::string param_file
int main()
Definition: adol-c_atom.cpp:46
#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 DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
PetscErrorCode DMMoFEMKSPSetComputeRHS(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set KSP right hand side evaluation function
Definition: DMMMoFEM.cpp:584
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:533
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMMoFEM.cpp:482
PetscErrorCode DMMoFEMKSPSetComputeOperators(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
Set KSP operators and push mofem finite element methods.
Definition: DMMMoFEM.cpp:625
#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.
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)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
static char help[]
Definition: poisson_prob.cpp:8
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.
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
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:313
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
boost::shared_ptr< VolEle > postFaceVolume
Poisson(moab::Core &mb_instance, MoFEM::Core &core, int order)
MoFEMErrorCode applyBc()
const int mpiRank
MoFEMErrorCode runAnalysis()
MatrixDouble invJac
std::string domainField
MPI_Comm mpiComm
boost::shared_ptr< VolEle > volPipelineLhs
MoFEM::Interface & mField
MoFEMErrorCode assembleSystem()
MoFEMErrorCode setupSystem()
MoFEMErrorCode setOutput()
MoFEMErrorCode setIntegrationRule()
Simple * simpleInterface
boost::shared_ptr< VolEle > volPipelineRhs
moab::Interface & moab
MoFEMErrorCode readMesh()
MoFEMErrorCode solveSystem()
Postprocess on face.
VolumeElementForcesAndSourcesCore VolEle