v0.13.1
poisson_2d_dis_galerkin.cpp
Go to the documentation of this file.
1/**
2 * \file poisson_2d_dis_galerkin.cpp
3 * \example poisson_2d_dis_galerkin.cpp
4 *
5 * Example of implementation for discontinuous Galerkin.
6 */
7
8#include <BasicFiniteElements.hpp>
9
10template <int DIM> struct ElementsAndOps {};
11
12template <> struct ElementsAndOps<2> {
17};
18
19constexpr int BASE_DIM = 1;
20constexpr int FIELD_DIM = 1;
21constexpr int SPACE_DIM = 2;
22
26
31
33
34static double penalty = 1e6;
35static double phi =
36 -1; // 1 - symmetric Nitsche, 0 - nonsymmetric, -1 antisymmetrica
37static double nitsche = 1;
38
40
45
46PetscBool is_test = PETSC_FALSE;
47
48auto u_exact = [](const double x, const double y, const double) {
49 if (is_test)
50 return x * x * y * y;
51 else
52 return cos(2 * x * M_PI) * cos(2 * y * M_PI);
53};
54
55auto u_grad_exact = [](const double x, const double y, const double) {
56 if (is_test)
57 return FTensor::Tensor1<double, 2>{2 * x * y * y, 2 * x * x * y};
58 else
60
61 -2 * M_PI * cos(2 * M_PI * y) * sin(2 * M_PI * x),
62 -2 * M_PI * cos(2 * M_PI * x) * sin(2 * M_PI * y)
63
64 };
65};
66
67auto source = [](const double x, const double y, const double) {
68 if (is_test)
69 return -(2 * x * x + 2 * y * y);
70 else
71 return 8 * M_PI * M_PI * cos(2 * x * M_PI) * cos(2 * y * M_PI);
72};
73
74using namespace MoFEM;
76
77static char help[] = "...\n\n";
78
80public:
82
83 // Declaration of the main function to run analysis
85
86private:
87 // Declaration of other main functions called in runProgram()
96
97 // MoFEM interfaces
100
101 // Field name and approximation order
102 std::string domainField;
103 int oRder;
104};
105
107 : domainField("U"), mField(m_field), oRder(4) {}
108
109//! [Read mesh]
112
115
116 // Only L2 field is set in this example. Two lines bellow forces simple
117 // interface to creat lower dimension (edge) elements, despite that fact that
118 // there is no field spanning on such elements. We need them for DG method.
121
123
125}
126//! [Read mesh]
127
128//! [Setup problem]
131
132 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &oRder, PETSC_NULL);
133 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-penalty", &penalty,
134 PETSC_NULL);
135 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-phi", &phi, PETSC_NULL);
136 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-nitsche", &nitsche,
137 PETSC_NULL);
138 PetscOptionsGetBool(PETSC_NULL, "", "-is_test", &is_test, PETSC_NULL);
139
140 MOFEM_LOG("WORLD", Sev::inform) << "Set order: " << oRder;
141 MOFEM_LOG("WORLD", Sev::inform) << "Set penalty: " << penalty;
142 MOFEM_LOG("WORLD", Sev::inform) << "Set phi: " << phi;
143 MOFEM_LOG("WORLD", Sev::inform) << "Set nitche: " << nitsche;
144 MOFEM_LOG("WORLD", Sev::inform) << "Set test: " << (is_test == PETSC_TRUE);
145
150
151 // This is only for debigging and experimentation, to see boundary and edge
152 // elements.
153 auto save_shared = [&](auto meshset, std::string prefix) {
155 auto file_name =
156 prefix + "_" +
157 boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk";
158 CHKERR mField.get_moab().write_file(file_name.c_str(), "VTK", "", &meshset,
159 1);
161 };
162
163 // CHKERR save_shared(simpleInterface->getBoundaryMeshSet(), "bdy");
164 // CHKERR save_shared(simpleInterface->getSkeletonMeshSet(), "skel");
165
167}
168//! [Setup problem]
169
170//! [Boundary condition]
174}
175
176//! [Boundary condition]
177
178//! [Assemble system]
181
182 auto pipeline_mng = mField.getInterface<PipelineManager>();
183
184 auto add_base_ops = [&](auto &pipeline) {
185 auto det_ptr = boost::make_shared<VectorDouble>();
186 auto jac_ptr = boost::make_shared<MatrixDouble>();
187 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
188 pipeline.push_back(new OpCalculateHOJacForFace(jac_ptr));
189 pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
190 pipeline.push_back(new OpSetInvJacL2ForFace(inv_jac_ptr));
191 };
192
193 add_base_ops(pipeline_mng->getOpDomainLhsPipeline());
194 pipeline_mng->getOpDomainLhsPipeline().push_back(new OpDomainGradGrad(
196 [](const double, const double, const double) { return 1; }));
197 pipeline_mng->getOpDomainRhsPipeline().push_back(
199
200 // Push operators to the Pipeline for Skeleton
201 auto side_fe_ptr = boost::make_shared<FaceSideEle>(mField);
202 add_base_ops(side_fe_ptr->getOpPtrVector());
203 side_fe_ptr->getOpPtrVector().push_back(
205
206 // Push operators to the Pipeline for Skeleton
207 pipeline_mng->getOpSkeletonLhsPipeline().push_back(
208 new OpL2LhsPenalty(side_fe_ptr));
209
210 // Push operators to the Pipeline for Boundary
211 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
212 new OpL2LhsPenalty(side_fe_ptr));
213 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
214 new OpL2BoundaryRhs(side_fe_ptr, u_exact));
215
217}
218//! [Assemble system]
219
220//! [Set integration rules]
223
224 auto rule_lhs = [](int, int, int p) -> int { return 2 * p; };
225 auto rule_rhs = [](int, int, int p) -> int { return 2 * p; };
226 auto rule_2 = [this](int, int, int) { return 2 * oRder; };
227
228 auto pipeline_mng = mField.getInterface<PipelineManager>();
229 CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule_lhs);
230 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule_rhs);
231
232 CHKERR pipeline_mng->setSkeletonLhsIntegrationRule(rule_2);
233 CHKERR pipeline_mng->setSkeletonRhsIntegrationRule(rule_2);
234 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(rule_2);
235 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(rule_2);
236
238}
239//! [Set integration rules]
240
241//! [Solve system]
244
245 auto pipeline_mng = mField.getInterface<PipelineManager>();
246
247 auto ksp_solver = pipeline_mng->createKSP();
248 CHKERR KSPSetFromOptions(ksp_solver);
249 CHKERR KSPSetUp(ksp_solver);
250
251 // Create RHS and solution vectors
252 auto dm = simpleInterface->getDM();
253 auto F = smartCreateDMVector(dm);
254 auto D = smartVectorDuplicate(F);
255
256 CHKERR KSPSolve(ksp_solver, F, D);
257
258 // Scatter result data on the mesh
259 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
260 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
261 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
262
264}
265//! [Solve system]
266
269
270 auto pipeline_mng = mField.getInterface<PipelineManager>();
271 pipeline_mng->getDomainRhsFE().reset();
272 pipeline_mng->getDomainLhsFE().reset();
273 pipeline_mng->getSkeletonRhsFE().reset();
274 pipeline_mng->getSkeletonLhsFE().reset();
275 pipeline_mng->getBoundaryRhsFE().reset();
276 pipeline_mng->getBoundaryLhsFE().reset();
277
278 auto rule = [](int, int, int p) -> int { return 2 * p; };
279 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
280 auto rule_2 = [this](int, int, int) { return 2 * oRder; };
281 CHKERR pipeline_mng->setSkeletonRhsIntegrationRule(rule_2);
282
283 auto add_base_ops = [&](auto &pipeline) {
284 auto det_ptr = boost::make_shared<VectorDouble>();
285 auto jac_ptr = boost::make_shared<MatrixDouble>();
286 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
287 pipeline.push_back(new OpCalculateHOJacForFace(jac_ptr));
288 pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
289 pipeline.push_back(new OpSetInvJacL2ForFace(inv_jac_ptr));
290 };
291
292 auto u_vals_ptr = boost::make_shared<VectorDouble>();
293 auto u_grad_ptr = boost::make_shared<MatrixDouble>();
294
295 add_base_ops(pipeline_mng->getOpDomainRhsPipeline());
296 pipeline_mng->getOpDomainRhsPipeline().push_back(
298 pipeline_mng->getOpDomainRhsPipeline().push_back(
300
301 enum NORMS { L2 = 0, SEMI_NORM, H1, LAST_NORM };
302 std::array<int, LAST_NORM> error_indices;
303 for (int i = 0; i != LAST_NORM; ++i)
304 error_indices[i] = i;
305 auto l2_vec = createSmartVectorMPI(
306 mField.get_comm(), (!mField.get_comm_rank()) ? LAST_NORM : 0, LAST_NORM);
307
308 auto error_op = new DomainEleOp(domainField, DomainEleOp::OPROW);
309 error_op->doWorkRhsHook = [&](DataOperator *op_ptr, int side, EntityType type,
312 auto o = static_cast<DomainEleOp *>(op_ptr);
313
315
316 if (const size_t nb_dofs = data.getIndices().size()) {
317
318 const int nb_integration_pts = o->getGaussPts().size2();
319 auto t_w = o->getFTensor0IntegrationWeight();
320 auto t_val = getFTensor0FromVec(*u_vals_ptr);
321 auto t_grad = getFTensor1FromMat<2>(*u_grad_ptr);
322 auto t_coords = o->getFTensor1CoordsAtGaussPts();
323
324 std::array<double, LAST_NORM> error;
325 std::fill(error.begin(), error.end(), 0);
326
327 for (int gg = 0; gg != nb_integration_pts; ++gg) {
328
329 const double alpha = t_w * o->getMeasure();
330 const double diff =
331 t_val - u_exact(t_coords(0), t_coords(1), t_coords(2));
332
333 auto t_exact_grad = u_grad_exact(t_coords(0), t_coords(1), t_coords(2));
334
335 const double diff_grad =
336 (t_grad(i) - t_exact_grad(i)) * (t_grad(i) - t_exact_grad(i));
337
338 error[L2] += alpha * pow(diff, 2);
339 error[SEMI_NORM] += alpha * diff_grad;
340
341 ++t_w;
342 ++t_val;
343 ++t_grad;
344 ++t_coords;
345 }
346
347 error[H1] = error[L2] + error[SEMI_NORM];
348
349 CHKERR VecSetValues(l2_vec, LAST_NORM, error_indices.data(), error.data(),
350 ADD_VALUES);
351 }
352
354 };
355
356 auto side_fe_ptr = boost::make_shared<FaceSideEle>(mField);
357 add_base_ops(side_fe_ptr->getOpPtrVector());
358 side_fe_ptr->getOpPtrVector().push_back(
360 std::array<VectorDouble, 2> side_vals;
361 std::array<double, 2> area_map;
362
363 auto side_vals_op = new DomainEleOp(domainField, DomainEleOp::OPROW);
364 side_vals_op->doWorkRhsHook = [&](DataOperator *op_ptr, int side,
368 auto o = static_cast<FaceSideOp *>(op_ptr);
369
370 const auto nb_in_loop = o->getFEMethod()->nInTheLoop;
371 area_map[nb_in_loop] = o->getMeasure();
372 side_vals[nb_in_loop] = *u_vals_ptr;
373 if (!nb_in_loop) {
374 area_map[1] = 0;
375 side_vals[1].clear();
376 }
377
379 };
380 side_fe_ptr->getOpPtrVector().push_back(side_vals_op);
381
382 auto do_work_rhs_error = [&](DataOperator *op_ptr, int side, EntityType type,
385 auto o = static_cast<BoundaryEleOp *>(op_ptr);
386
387 CHKERR o->loopSideFaces("dFE", *side_fe_ptr);
388 const auto in_the_loop = side_fe_ptr->nInTheLoop;
389
390#ifndef NDEBUG
391 const std::array<std::string, 2> ele_type_name = {"BOUNDARY", "SKELETON"};
392 MOFEM_LOG("SELF", Sev::noisy)
393 << "do_work_rhs_error in_the_loop " << ele_type_name[in_the_loop];
394#endif
395
396 const double s = o->getMeasure() / (area_map[0] + area_map[1]);
397 const double p = penalty * s;
398
399 constexpr std::array<int, 2> sign_array{1, -1};
400
401 std::array<double, LAST_NORM> error;
402 std::fill(error.begin(), error.end(), 0);
403
404 const int nb_integration_pts = o->getGaussPts().size2();
405
406 if (!in_the_loop) {
407 side_vals[1].resize(nb_integration_pts, false);
408 auto t_coords = o->getFTensor1CoordsAtGaussPts();
409 auto t_val_m = getFTensor0FromVec(side_vals[1]);
410 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
411 t_val_m = u_exact(t_coords(0), t_coords(1), t_coords(2));
412 ++t_coords;
413 ++t_val_m;
414 }
415 }
416
417 auto t_val_p = getFTensor0FromVec(side_vals[0]);
418 auto t_val_m = getFTensor0FromVec(side_vals[1]);
419 auto t_w = o->getFTensor0IntegrationWeight();
420
421 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
422
423 const double alpha = o->getMeasure() * t_w;
424 const auto diff = t_val_p - t_val_m;
425 error[SEMI_NORM] += alpha * p * diff * diff;
426
427 ++t_w;
428 ++t_val_p;
429 ++t_val_m;
430 }
431
432
433 error[H1] = error[SEMI_NORM];
434 CHKERR VecSetValues(l2_vec, LAST_NORM, error_indices.data(), error.data(),
435 ADD_VALUES);
436
438 };
439
440 auto skeleton_error_op = new BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE);
441 skeleton_error_op->doWorkRhsHook = do_work_rhs_error;
442 auto boundary_error_op = new BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE);
443 boundary_error_op->doWorkRhsHook = do_work_rhs_error;
444
445 pipeline_mng->getOpDomainRhsPipeline().push_back(error_op);
446 pipeline_mng->getOpSkeletonLhsPipeline().push_back(skeleton_error_op);
447 pipeline_mng->getOpBoundaryRhsPipeline().push_back(boundary_error_op);
448
449 CHKERR pipeline_mng->loopFiniteElements();
450
451 CHKERR VecAssemblyBegin(l2_vec);
452 CHKERR VecAssemblyEnd(l2_vec);
453
454 if (mField.get_comm_rank() == 0) {
455 const double *array;
456 CHKERR VecGetArrayRead(l2_vec, &array);
457 MOFEM_LOG_C("SELF", Sev::inform, "Error Norm L2 %6.4e",
458 std::sqrt(array[L2]));
459 MOFEM_LOG_C("SELF", Sev::inform, "Error Norm Energetic %6.4e",
460 std::sqrt(array[SEMI_NORM]));
461 MOFEM_LOG_C("SELF", Sev::inform, "Error Norm H1 %6.4e",
462 std::sqrt(array[H1]));
463
464 if(is_test) {
465 constexpr double eps = 1e-12;
466 if (std::sqrt(array[H1]) > eps)
467 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID, "Error is too big");
468 }
469
470 CHKERR VecRestoreArrayRead(l2_vec, &array);
471 const MoFEM::Problem *problem_ptr;
473 MOFEM_LOG_C("SELF", Sev::inform, "Nb. DOFs %d",
474 problem_ptr->getNbDofsRow());
475 }
476
477
478
480}
481
482//! [Output results]
485
486 auto pipeline_mng = mField.getInterface<PipelineManager>();
487 pipeline_mng->getDomainLhsFE().reset();
488 pipeline_mng->getSkeletonRhsFE().reset();
489 pipeline_mng->getSkeletonLhsFE().reset();
490 pipeline_mng->getBoundaryRhsFE().reset();
491 pipeline_mng->getBoundaryLhsFE().reset();
492
493 auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
494
495 auto u_ptr = boost::make_shared<VectorDouble>();
496 post_proc_fe->getOpPtrVector().push_back(
498
500
501 post_proc_fe->getOpPtrVector().push_back(
502
503 new OpPPMap(
504
505 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
506
507 {{"U", u_ptr}},
508
509 {},
510
511 {},
512
513 {})
514
515 );
516
517 pipeline_mng->getDomainRhsFE() = post_proc_fe;
518 CHKERR pipeline_mng->loopFiniteElements();
519 CHKERR post_proc_fe->writeFile("out_result.h5m");
520
522}
523//! [Output results]
524
525//! [Run program]
528
537
539}
540//! [Run program]
541
542//! [Main]
543int main(int argc, char *argv[]) {
544
545 // Initialisation of MoFEM/PETSc and MOAB data structures
546 const char param_file[] = "param_file.petsc";
548
549 // Error handling
550 try {
551 // Register MoFEM discrete manager in PETSc
552 DMType dm_name = "DMMOFEM";
553 CHKERR DMRegister_MoFEM(dm_name);
554
555 // Create MOAB instance
556 moab::Core mb_instance; // mesh database
557 moab::Interface &moab = mb_instance; // mesh database interface
558
559 // Create MoFEM instance
560 MoFEM::Core core(moab); // finite element database
561 MoFEM::Interface &m_field = core; // finite element interface
562
563 // Run the main analysis
564 Poisson2DiscontGalerkin poisson_problem(m_field);
565 CHKERR poisson_problem.runProgram();
566 }
568
569 // Finish work: cleaning memory, getting statistics, etc.
571
572 return 0;
573}
574//! [Main]
static Index< 'o', 3 > o
static Index< 'p', 3 > p
std::string param_file
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:304
MoFEM::FaceElementForcesAndSourcesCore FaceEle
static const double eps
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
@ NOSPACE
Definition: definitions.h:83
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
#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
ElementsAndOps< SPACE_DIM >::FaceSideEle FaceSideEle
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:470
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:373
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
FTensor::Index< 'i', SPACE_DIM > i
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
Definition: helmholtz.cpp:27
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
auto createSmartVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
CoreTmp< 0 > Core
Definition: Core.hpp:1086
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
const double D
diffusivity
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
int main(int argc, char *argv[])
[Run program]
EntitiesFieldData::EntData EntData
auto u_grad_exact
static char help[]
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::BoundaryEleOp BoundaryEleOp
static double nitsche
constexpr int FIELD_DIM
constexpr int BASE_DIM
DomainEle::UserDataOperator DomainEleOp
static double penalty
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< BASE_DIM, FIELD_DIM, SPACE_DIM > OpDomainGradGrad
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, FIELD_DIM > OpDomainSource
static double phi
PetscBool is_test
virtual moab::Interface & get_moab()=0
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
base operator to do operations at Gauss Pt. level
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
Base face element used to integrate on skeleton.
@ OPROW
operator doWork function is executed on FE rows
@ OPSPACE
operator do Work is execute on space data
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.
Definition: PostProc.hpp:166
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Definition: PostProc.hpp:120
keeps basic data about problem
DofIdx getNbDofsRow() const
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
bool & getAddBoundaryFE()
Get the addSkeletonFE.
Definition: Simple.hpp:408
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:374
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:290
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:806
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:304
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:589
bool & getAddSkeletonFE()
Get the addSkeletonFE.
Definition: Simple.hpp:397
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:749
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MoFEMErrorCode boundaryCondition()
[Setup problem]
MoFEMErrorCode checkResults()
[Solve system]
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode setIntegrationRules()
[Assemble system]
MoFEMErrorCode assembleSystem()
[Boundary condition]
Poisson2DiscontGalerkin(MoFEM::Interface &m_field)
MoFEMErrorCode solveSystem()
[Set integration rules]
MoFEMErrorCode readMesh()
[Read mesh]
MoFEMErrorCode outputResults()
[Output results]
MoFEMErrorCode runProgram()
[Output results]
Operator tp collect data from elements on the side of Edge/Face.
Opator tp evaluate Dirichlet boundary conditions using DG.