v0.14.0
Loading...
Searching...
No Matches
standard_poisson.cpp
Go to the documentation of this file.
1#include <stdlib.h>
3
4using namespace MoFEM;
5// using namespace Poisson2DNonhomogeneousOperators;
6
7inline double sqr(double x) { return x * x; }
8
9constexpr int SPACE_DIM = 2;
12
15
18
21
22static char help[] = "...\n\n";
23
25public:
27
28 // Declaration of the main function to run analysis
30
31private:
32 // Declaration of other main functions called in runProgram()
42
43 //! [Analytical function]
44 static double analyticalFunction(const double x, const double y,
45 const double z) {
46 return exp(-100. * (sqr(x) + sqr(y))) * cos(M_PI * x) * cos(M_PI * y);
47 }
48 //! [Analytical function]
49
50 //! [Analytical function gradient]
51 static VectorDouble analyticalFunctionGrad(const double x, const double y,
52 const double z) {
53 VectorDouble res;
54 res.resize(2);
55 res[0] = -exp(-100. * (sqr(x) + sqr(y))) *
56 (200. * x * cos(M_PI * x) + M_PI * sin(M_PI * x)) * cos(M_PI * y);
57 res[1] = -exp(-100. * (sqr(x) + sqr(y))) *
58 (200. * y * cos(M_PI * y) + M_PI * sin(M_PI * y)) * cos(M_PI * x);
59 return res;
60 }
61 //! [Analytical function gradient]
62
63 //! [Source function]
64 static double sourceFunction(const double x, const double y, const double z) {
65 return -exp(-100. * (sqr(x) + sqr(y))) *
66 (400. * M_PI *
67 (x * cos(M_PI * y) * sin(M_PI * x) +
68 y * cos(M_PI * x) * sin(M_PI * y)) +
69 2. * (20000. * (sqr(x) + sqr(y)) - 200. - sqr(M_PI)) *
70 cos(M_PI * x) * cos(M_PI * y));
71 }
72 //! [Source function]
73
74 struct CommonData {
75 boost::shared_ptr<VectorDouble> approxVals;
76 boost::shared_ptr<MatrixDouble> approxValsGrad;
78
80 };
81
82 boost::shared_ptr<CommonData> commonDataPtr;
83
84 struct OpError : public DomainEleOp {
85 std::string domainField;
86 boost::shared_ptr<CommonData> commonDataPtr;
88 OpError(std::string domain_field,
89 boost::shared_ptr<CommonData> &common_data_ptr,
90 MoFEM::Interface &m_field)
91 : DomainEleOp(domain_field, OPROW), domainField(domain_field),
92 commonDataPtr(common_data_ptr), mField(m_field) {
93 std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
94 doEntities[MBTRI] = doEntities[MBQUAD] = true;
95 }
96 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
97 };
98
99 // Main interfaces
102
103 // Field name and approximation order
104 std::string domainField;
105 int oRder;
106
107 // Object to mark boundary entities for the assembling of domain elements
108 boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
109
110 // Boundary entities marked for fieldsplit (block) solver - optional
112};
113
115 : domainField("U"), mField(m_field) {}
116
119
123
125}
126
129
134
135 int oRder = 2;
136 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &oRder, PETSC_NULL);
138
140
142}
143
146
147 auto bc_mng = mField.getInterface<BcManager>();
148
149 // Remove BCs from blockset name "BOUNDARY_CONDITION" or SETU, note that you
150 // can use regular expression to put list of blocksets;
151 CHKERR bc_mng->removeBlockDOFsOnEntities<BcScalarMeshsetType<BLOCKSET>>(
152 simpleInterface->getProblemName(), "BOUNDARY", std::string(domainField),
153 true);
154
156}
157
158//! [Create common data]
161 commonDataPtr = boost::make_shared<CommonData>();
162 PetscInt ghosts[2] = {0, 1};
163 if (!mField.get_comm_rank())
164 commonDataPtr->petscVec =
165 createGhostVector(mField.get_comm(), 2, 2, 0, ghosts);
166 else
167 commonDataPtr->petscVec =
168 createGhostVector(mField.get_comm(), 0, 2, 2, ghosts);
169 commonDataPtr->approxVals = boost::make_shared<VectorDouble>();
170 commonDataPtr->approxValsGrad = boost::make_shared<MatrixDouble>();
172}
173//! [Create common data]
174
177
178 auto pipeline_mng = mField.getInterface<PipelineManager>();
179
180 { // Push operators to the Pipeline that is responsible for calculating LHS of
181 // domain elements
182 CHKERR AddHOOps<2, 2, 2>::add(pipeline_mng->getOpDomainLhsPipeline(), {H1});
183 auto unity = [](const double, const double, const double) { return 1; };
184 pipeline_mng->getOpDomainLhsPipeline().push_back(
186 }
187
188 { // Push operators to the Pipeline that is responsible for calculating RHS of
189 // domain elements
190 pipeline_mng->getOpDomainRhsPipeline().push_back(
192 }
193
195}
196
199
200 auto pipeline_mng = mField.getInterface<PipelineManager>();
201
202 auto domain_rule_lhs = [](int, int, int p) -> int { return 2 * (p - 1); };
203 auto domain_rule_rhs = [](int, int, int p) -> int { return 2 * (p - 1); };
204 CHKERR pipeline_mng->setDomainLhsIntegrationRule(domain_rule_lhs);
205 CHKERR pipeline_mng->setDomainRhsIntegrationRule(domain_rule_rhs);
206
207 auto boundary_rule_lhs = [](int, int, int p) -> int { return 2 * p; };
208 auto boundary_rule_rhs = [](int, int, int p) -> int { return 2 * p; };
209 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(boundary_rule_lhs);
210 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(boundary_rule_rhs);
211
213}
214
217
218 auto pipeline_mng = mField.getInterface<PipelineManager>();
219
220 auto ksp_solver = pipeline_mng->createKSP();
221 CHKERR KSPSetFromOptions(ksp_solver);
222
223 // Create RHS and solution vectors
224 auto dm = simpleInterface->getDM();
225 auto F = createDMVector(dm);
226 auto D = vectorDuplicate(F);
227
228 CHKERR KSPSetUp(ksp_solver);
229
230 // Solve the system
231 CHKERR KSPSolve(ksp_solver, F, D);
232
233 // Scatter result data on the mesh
234 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
235 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
236 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
237
239}
240
241//! [Check error]
245 pipeline_mng->getDomainLhsFE().reset();
246 pipeline_mng->getDomainRhsFE().reset();
247 pipeline_mng->getOpDomainRhsPipeline().clear();
248
250
251 pipeline_mng->getOpDomainRhsPipeline().push_back(
253 pipeline_mng->getOpDomainRhsPipeline().push_back(
255 commonDataPtr->approxValsGrad));
256 pipeline_mng->getOpDomainRhsPipeline().push_back(
258
259 CHKERR VecZeroEntries(commonDataPtr->petscVec);
260
261 CHKERR pipeline_mng->loopFiniteElements();
262
263 CHKERR VecAssemblyBegin(commonDataPtr->petscVec);
264 CHKERR VecAssemblyEnd(commonDataPtr->petscVec);
265 const double *array;
266 CHKERR VecGetArrayRead(commonDataPtr->petscVec, &array);
267 MOFEM_LOG("EXAMPLE", Sev::inform)
268 << "Global error L2 norm: " << std::sqrt(array[0]);
269 MOFEM_LOG("EXAMPLE", Sev::inform)
270 << "Global error H1 seminorm: " << std::sqrt(array[1]);
271 CHKERR VecRestoreArrayRead(commonDataPtr->petscVec, &array);
272
274}
275//! [Check error]
276
279
280 auto pipeline_mng = mField.getInterface<PipelineManager>();
281 pipeline_mng->getDomainLhsFE().reset();
282 pipeline_mng->getBoundaryLhsFE().reset();
283 pipeline_mng->getDomainRhsFE().reset();
284 pipeline_mng->getBoundaryRhsFE().reset();
285
287
288 auto post_proc_fe = boost::make_shared<PostProcFaceEle>(mField);
289
290 CHKERR AddHOOps<2, 2, 2>::add(post_proc_fe->getOpPtrVector(), {H1});
291
292 post_proc_fe->getOpPtrVector().push_back(
294 post_proc_fe->getOpPtrVector().push_back(
296 commonDataPtr->approxValsGrad));
297
298 post_proc_fe->getOpPtrVector().push_back(new OpPPMap(
299 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
300 {{domainField, commonDataPtr->approxVals}},
301 {{domainField + "_GRAD", commonDataPtr->approxValsGrad}}, {}, {}));
302 pipeline_mng->getDomainRhsFE() = post_proc_fe;
303
304 CHKERR pipeline_mng->loopFiniteElements();
305 CHKERR post_proc_fe->writeFile("out_result.h5m");
306
308}
309
312
322
324}
325
326int main(int argc, char *argv[]) {
327
328 // Initialisation of MoFEM/PETSc and MOAB data structures
329 const char param_file[] = "param_file.petsc";
331
332 // Add logging channel for example problem
333 auto core_log = logging::core::get();
334 core_log->add_sink(
336 LogManager::setLog("EXAMPLE");
337 MOFEM_LOG_TAG("EXAMPLE", "StandardPoisson");
338
339 // Error handling
340 try {
341 // Register MoFEM discrete manager in PETSc
342 DMType dm_name = "DMMOFEM";
343 CHKERR DMRegister_MoFEM(dm_name);
344
345 // Create MOAB instance
346 moab::Core mb_instance; // mesh database
347 moab::Interface &moab = mb_instance; // mesh database interface
348
349 // Create MoFEM instance
350 MoFEM::Core core(moab); // finite element database
351 MoFEM::Interface &m_field = core; // finite element interface
352
353 // Run the main analysis
354 StandardPoisson poisson_problem(m_field);
355 CHKERR poisson_problem.runProgram();
356 }
358
359 // Finish work: cleaning memory, getting statistics, etc.
361
362 return 0;
363}
364
365//! [OpError]
367 EntData &data) {
369 const int nb_integration_pts = getGaussPts().size2();
370 const double area = getMeasure();
371 auto t_w = getFTensor0IntegrationWeight();
372 auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
373 auto t_val_grad = getFTensor1FromMat<2>(*(commonDataPtr->approxValsGrad));
374 auto t_coords = getFTensor1CoordsAtGaussPts();
377
378 double error_l2 = 0;
379 double error_h1 = 0;
380
381 for (int gg = 0; gg != nb_integration_pts; ++gg) {
382 const double alpha = t_w * area;
383
384 double diff = t_val - StandardPoisson::analyticalFunction(
385 t_coords(0), t_coords(1), t_coords(2));
386 error_l2 += alpha * sqr(diff);
387
389 t_coords(0), t_coords(1), t_coords(2));
390 auto t_fun_grad = getFTensor1FromArray<2, 2>(vec);
391 t_diff(i) = t_val_grad(i) - t_fun_grad(i);
392
393 error_h1 += alpha * t_diff(i) * t_diff(i);
394
395 ++t_w;
396 ++t_val;
397 ++t_val_grad;
398 ++t_coords;
399 }
400
401 int index = CommonData::ERROR_L2_NORM;
402 constexpr std::array<int, 2> indices = {CommonData::ERROR_L2_NORM,
404 std::array<double, 2> values;
405 values[0] = error_l2;
406 values[1] = error_h1;
407 CHKERR VecSetValues(commonDataPtr->petscVec, 2, indices.data(), values.data(),
408 ADD_VALUES);
410}
411//! [OpError]
static Index< 'p', 3 > p
std::string param_file
int main()
Definition: adol-c_atom.cpp:46
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
#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
#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
@ 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 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
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
FTensor::Index< 'i', SPACE_DIM > i
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
Definition: helmholtz.cpp:25
double D
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)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
auto createGhostVector(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[])
Create smart ghost vector.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
double sqr(double x)
static char help[]
constexpr int SPACE_DIM
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpDomainSource
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =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
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
@ OPROW
operator doWork function is executed on FE rows
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
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.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEM::FaceElementForcesAndSourcesCore FaceEle
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
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
boost::shared_ptr< VectorDouble > approxVals
SmartPetscObj< Vec > petscVec
boost::shared_ptr< MatrixDouble > approxValsGrad
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[OpError]
MoFEM::Interface & mField
OpError(std::string domain_field, boost::shared_ptr< CommonData > &common_data_ptr, MoFEM::Interface &m_field)
MoFEMErrorCode outputResults()
[Check error]
MoFEMErrorCode setIntegrationRules()
static double analyticalFunction(const double x, const double y, const double z)
[Analytical function]
boost::shared_ptr< CommonData > commonDataPtr
StandardPoisson(MoFEM::Interface &m_field)
static double sourceFunction(const double x, const double y, const double z)
[Analytical function gradient]
MoFEMErrorCode createCommonData()
[Create common data]
MoFEMErrorCode checkError()
[Check error]
std::string domainField
MoFEMErrorCode runProgram()
MoFEMErrorCode solveSystem()
MoFEMErrorCode assembleSystem()
[Create common data]
Range boundaryEntitiesForFieldsplit
MoFEMErrorCode setupProblem()
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode readMesh()
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
MoFEM::Interface & mField
static VectorDouble analyticalFunctionGrad(const double x, const double y, const double z)
[Analytical function]