v0.13.0
poisson_prob.cpp
Go to the documentation of this file.
1 #include <stdlib.h>
2 #include <BasicFiniteElements.hpp>
3 #include <PoissonOperators.hpp>
4 
5 using namespace MoFEM;
6 using namespace PoissonOps;
7 
8 static char help[] = "...\n\n";
9 
10 struct Poisson {
11 public:
12  Poisson(moab::Core &mb_instance, MoFEM::Core &core, int order);
13 
14  MoFEMErrorCode runAnalysis();
15 
16 private:
17  MoFEMErrorCode readMesh();
18  MoFEMErrorCode setupSystem();
19  MoFEMErrorCode setIntegrationRule();
20  MoFEMErrorCode assembleSystem();
21  MoFEMErrorCode applyBc();
22  MoFEMErrorCode solveSystem();
23  MoFEMErrorCode setOutput();
24 
25  boost::shared_ptr<VolEle> volPipelineLhs;
26  boost::shared_ptr<VolEle> volPipelineRhs;
27  boost::shared_ptr<VolEle> postFaceVolume;
28 
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 
46 Poisson::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 
197  CHKERR readMesh();
199  CHKERR applyBc();
200 
202 
203  CHKERR setOutput();
204 
206 
208 }
209 
210 int main(int argc, char *argv[]) {
211  const char param_file[] = "param_file.petsc";
212  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
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  }
227  CATCH_ERRORS;
229  return 0;
230 }
static Index< 'p', 3 > p
std::string param_file
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ H1
continuous field
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ BLOCKSET
Definition: definitions.h:161
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
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:595
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:544
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMMoFEM.cpp:493
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:636
#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.
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
CoreTmp< 0 > Core
Definition: Core.hpp:1096
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
int main(int argc, char *argv[])
static char help[]
Definition: poisson_prob.cpp:8
virtual moab::Interface & get_moab()=0
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:33
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:293
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:216
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:715
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:230
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:508
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:668
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:340
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:319
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