v0.14.0
poisson_2d_dis_galerkin.cpp

Example of implementation for discontinuous Galerkin.

/**
* \file poisson_2d_dis_galerkin.cpp
* \example poisson_2d_dis_galerkin.cpp
*
* Example of implementation for discontinuous Galerkin.
*/
constexpr int BASE_DIM = 1;
constexpr int FIELD_DIM = 1;
constexpr int SPACE_DIM = 2;
using BoundaryEle =
using FaceSideEle =
using PostProcEle = PostProcBrokenMeshInMoab<DomainEle>;
static double penalty = 1e6;
static double phi =
-1; // 1 - symmetric Nitsche, 0 - nonsymmetric, -1 antisymmetrica
static double nitsche = 1;
using OpDomainGradGrad = FormsIntegrators<DomainEleOp>::Assembly<
using OpDomainSource = FormsIntegrators<DomainEleOp>::Assembly<
PETSC>::LinearForm<GAUSS>::OpSource<BASE_DIM, FIELD_DIM>;
PetscBool is_test = PETSC_FALSE;
auto u_exact = [](const double x, const double y, const double) {
if (is_test)
return x * x * y * y;
else
return cos(2 * x * M_PI) * cos(2 * y * M_PI);
};
auto u_grad_exact = [](const double x, const double y, const double) {
if (is_test)
return FTensor::Tensor1<double, 2>{2 * x * y * y, 2 * x * x * y};
else
-2 * M_PI * cos(2 * M_PI * y) * sin(2 * M_PI * x),
-2 * M_PI * cos(2 * M_PI * x) * sin(2 * M_PI * y)
};
};
auto source = [](const double x, const double y, const double) {
if (is_test)
return -(2 * x * x + 2 * y * y);
else
return 8 * M_PI * M_PI * cos(2 * x * M_PI) * cos(2 * y * M_PI);
};
using namespace MoFEM;
static char help[] = "...\n\n";
public:
// Declaration of the main function to run analysis
MoFEMErrorCode runProgram();
private:
// Declaration of other main functions called in runProgram()
MoFEMErrorCode readMesh();
MoFEMErrorCode setupProblem();
MoFEMErrorCode boundaryCondition();
MoFEMErrorCode assembleSystem();
MoFEMErrorCode setIntegrationRules();
MoFEMErrorCode solveSystem();
MoFEMErrorCode checkResults();
MoFEMErrorCode outputResults();
// MoFEM interfaces
Simple *simpleInterface;
// Field name and approximation order
std::string domainField;
int oRder;
};
: domainField("U"), mField(m_field), oRder(4) {}
//! [Read mesh]
// Only L2 field is set in this example. Two lines bellow forces simple
// interface to creat lower dimension (edge) elements, despite that fact that
// there is no field spanning on such elements. We need them for DG method.
}
//! [Read mesh]
//! [Setup problem]
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &oRder, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-penalty", &penalty,
PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-phi", &phi, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-nitsche", &nitsche,
PETSC_NULL);
PetscOptionsGetBool(PETSC_NULL, "", "-is_test", &is_test, PETSC_NULL);
MOFEM_LOG("WORLD", Sev::inform) << "Set order: " << oRder;
MOFEM_LOG("WORLD", Sev::inform) << "Set penalty: " << penalty;
MOFEM_LOG("WORLD", Sev::inform) << "Set phi: " << phi;
MOFEM_LOG("WORLD", Sev::inform) << "Set nitche: " << nitsche;
MOFEM_LOG("WORLD", Sev::inform) << "Set test: " << (is_test == PETSC_TRUE);
// This is only for debigging and experimentation, to see boundary and edge
// elements.
auto save_shared = [&](auto meshset, std::string prefix) {
auto file_name =
prefix + "_" +
boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk";
CHKERR mField.get_moab().write_file(file_name.c_str(), "VTK", "", &meshset,
1);
};
// CHKERR save_shared(simpleInterface->getBoundaryMeshSet(), "bdy");
// CHKERR save_shared(simpleInterface->getSkeletonMeshSet(), "skel");
}
//! [Setup problem]
//! [Boundary condition]
}
//! [Boundary condition]
//! [Assemble system]
auto pipeline_mng = mField.getInterface<PipelineManager>();
auto add_base_ops = [&](auto &pipeline) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(new OpCalculateHOJacForFace(jac_ptr));
pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(new OpSetInvJacL2ForFace(inv_jac_ptr));
};
add_base_ops(pipeline_mng->getOpDomainLhsPipeline());
pipeline_mng->getOpDomainLhsPipeline().push_back(new OpDomainGradGrad(
[](const double, const double, const double) { return 1; }));
pipeline_mng->getOpDomainRhsPipeline().push_back(
// Push operators to the Pipeline for Skeleton
auto side_fe_ptr = boost::make_shared<FaceSideEle>(mField);
add_base_ops(side_fe_ptr->getOpPtrVector());
side_fe_ptr->getOpPtrVector().push_back(
// Push operators to the Pipeline for Skeleton
pipeline_mng->getOpSkeletonLhsPipeline().push_back(
new OpL2LhsPenalty(side_fe_ptr));
// Push operators to the Pipeline for Boundary
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
new OpL2LhsPenalty(side_fe_ptr));
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
new OpL2BoundaryRhs(side_fe_ptr, u_exact));
}
//! [Assemble system]
//! [Set integration rules]
auto rule_lhs = [](int, int, int p) -> int { return 2 * p; };
auto rule_rhs = [](int, int, int p) -> int { return 2 * p; };
auto rule_2 = [this](int, int, int) { return 2 * oRder; };
auto pipeline_mng = mField.getInterface<PipelineManager>();
CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule_lhs);
CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule_rhs);
CHKERR pipeline_mng->setSkeletonLhsIntegrationRule(rule_2);
CHKERR pipeline_mng->setSkeletonRhsIntegrationRule(rule_2);
CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(rule_2);
CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(rule_2);
}
//! [Set integration rules]
//! [Solve system]
auto pipeline_mng = mField.getInterface<PipelineManager>();
auto ksp_solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(ksp_solver);
CHKERR KSPSetUp(ksp_solver);
// Create RHS and solution vectors
auto dm = simpleInterface->getDM();
auto F = createDMVector(dm);
auto D = vectorDuplicate(F);
CHKERR KSPSolve(ksp_solver, F, D);
// Scatter result data on the mesh
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
}
//! [Solve system]
auto pipeline_mng = mField.getInterface<PipelineManager>();
pipeline_mng->getDomainRhsFE().reset();
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getSkeletonRhsFE().reset();
pipeline_mng->getSkeletonLhsFE().reset();
pipeline_mng->getBoundaryRhsFE().reset();
pipeline_mng->getBoundaryLhsFE().reset();
auto rule = [](int, int, int p) -> int { return 2 * p; };
CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
auto rule_2 = [this](int, int, int) { return 2 * oRder; };
CHKERR pipeline_mng->setSkeletonRhsIntegrationRule(rule_2);
auto add_base_ops = [&](auto &pipeline) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(new OpCalculateHOJacForFace(jac_ptr));
pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(new OpSetInvJacL2ForFace(inv_jac_ptr));
};
auto u_vals_ptr = boost::make_shared<VectorDouble>();
auto u_grad_ptr = boost::make_shared<MatrixDouble>();
add_base_ops(pipeline_mng->getOpDomainRhsPipeline());
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
enum NORMS { L2 = 0, SEMI_NORM, H1, LAST_NORM };
std::array<int, LAST_NORM> error_indices;
for (int i = 0; i != LAST_NORM; ++i)
error_indices[i] = i;
auto l2_vec = createVectorMPI(
mField.get_comm(), (!mField.get_comm_rank()) ? LAST_NORM : 0, LAST_NORM);
auto error_op = new DomainEleOp(domainField, DomainEleOp::OPROW);
error_op->doWorkRhsHook = [&](DataOperator *op_ptr, int side, EntityType type,
auto o = static_cast<DomainEleOp *>(op_ptr);
if (const size_t nb_dofs = data.getIndices().size()) {
const int nb_integration_pts = o->getGaussPts().size2();
auto t_w = o->getFTensor0IntegrationWeight();
auto t_val = getFTensor0FromVec(*u_vals_ptr);
auto t_grad = getFTensor1FromMat<2>(*u_grad_ptr);
auto t_coords = o->getFTensor1CoordsAtGaussPts();
std::array<double, LAST_NORM> error;
std::fill(error.begin(), error.end(), 0);
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = t_w * o->getMeasure();
const double diff =
t_val - u_exact(t_coords(0), t_coords(1), t_coords(2));
auto t_exact_grad = u_grad_exact(t_coords(0), t_coords(1), t_coords(2));
const double diff_grad =
(t_grad(i) - t_exact_grad(i)) * (t_grad(i) - t_exact_grad(i));
error[L2] += alpha * pow(diff, 2);
error[SEMI_NORM] += alpha * diff_grad;
++t_w;
++t_val;
++t_grad;
++t_coords;
}
error[H1] = error[L2] + error[SEMI_NORM];
CHKERR VecSetValues(l2_vec, LAST_NORM, error_indices.data(), error.data(),
ADD_VALUES);
}
};
auto side_fe_ptr = boost::make_shared<FaceSideEle>(mField);
add_base_ops(side_fe_ptr->getOpPtrVector());
side_fe_ptr->getOpPtrVector().push_back(
std::array<VectorDouble, 2> side_vals;
std::array<double, 2> area_map;
auto side_vals_op = new DomainEleOp(domainField, DomainEleOp::OPROW);
side_vals_op->doWorkRhsHook = [&](DataOperator *op_ptr, int side,
EntityType type,
auto o = static_cast<FaceSideOp *>(op_ptr);
const auto nb_in_loop = o->getFEMethod()->nInTheLoop;
area_map[nb_in_loop] = o->getMeasure();
side_vals[nb_in_loop] = *u_vals_ptr;
if (!nb_in_loop) {
area_map[1] = 0;
side_vals[1].clear();
}
};
side_fe_ptr->getOpPtrVector().push_back(side_vals_op);
auto do_work_rhs_error = [&](DataOperator *op_ptr, int side, EntityType type,
auto o = static_cast<BoundaryEleOp *>(op_ptr);
CHKERR o->loopSideFaces("dFE", *side_fe_ptr);
const auto in_the_loop = side_fe_ptr->nInTheLoop;
#ifndef NDEBUG
const std::array<std::string, 2> ele_type_name = {"BOUNDARY", "SKELETON"};
MOFEM_LOG("SELF", Sev::noisy)
<< "do_work_rhs_error in_the_loop " << ele_type_name[in_the_loop];
#endif
const double s = o->getMeasure() / (area_map[0] + area_map[1]);
const double p = penalty * s;
constexpr std::array<int, 2> sign_array{1, -1};
std::array<double, LAST_NORM> error;
std::fill(error.begin(), error.end(), 0);
const int nb_integration_pts = o->getGaussPts().size2();
if (!in_the_loop) {
side_vals[1].resize(nb_integration_pts, false);
auto t_coords = o->getFTensor1CoordsAtGaussPts();
auto t_val_m = getFTensor0FromVec(side_vals[1]);
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
t_val_m = u_exact(t_coords(0), t_coords(1), t_coords(2));
++t_coords;
++t_val_m;
}
}
auto t_val_p = getFTensor0FromVec(side_vals[0]);
auto t_val_m = getFTensor0FromVec(side_vals[1]);
auto t_w = o->getFTensor0IntegrationWeight();
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = o->getMeasure() * t_w;
const auto diff = t_val_p - t_val_m;
error[SEMI_NORM] += alpha * p * diff * diff;
++t_w;
++t_val_p;
++t_val_m;
}
error[H1] = error[SEMI_NORM];
CHKERR VecSetValues(l2_vec, LAST_NORM, error_indices.data(), error.data(),
ADD_VALUES);
};
auto skeleton_error_op = new BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE);
skeleton_error_op->doWorkRhsHook = do_work_rhs_error;
auto boundary_error_op = new BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE);
boundary_error_op->doWorkRhsHook = do_work_rhs_error;
pipeline_mng->getOpDomainRhsPipeline().push_back(error_op);
pipeline_mng->getOpSkeletonLhsPipeline().push_back(skeleton_error_op);
pipeline_mng->getOpBoundaryRhsPipeline().push_back(boundary_error_op);
CHKERR pipeline_mng->loopFiniteElements();
CHKERR VecAssemblyBegin(l2_vec);
CHKERR VecAssemblyEnd(l2_vec);
if (mField.get_comm_rank() == 0) {
const double *array;
CHKERR VecGetArrayRead(l2_vec, &array);
MOFEM_LOG_C("SELF", Sev::inform, "Error Norm L2 %6.4e",
std::sqrt(array[L2]));
MOFEM_LOG_C("SELF", Sev::inform, "Error Norm Energetic %6.4e",
std::sqrt(array[SEMI_NORM]));
MOFEM_LOG_C("SELF", Sev::inform, "Error Norm H1 %6.4e",
std::sqrt(array[H1]));
if(is_test) {
constexpr double eps = 1e-12;
if (std::sqrt(array[H1]) > eps)
SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID, "Error is too big");
}
CHKERR VecRestoreArrayRead(l2_vec, &array);
const MoFEM::Problem *problem_ptr;
MOFEM_LOG_C("SELF", Sev::inform, "Nb. DOFs %d",
problem_ptr->getNbDofsRow());
}
}
//! [Output results]
auto pipeline_mng = mField.getInterface<PipelineManager>();
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getSkeletonRhsFE().reset();
pipeline_mng->getSkeletonLhsFE().reset();
pipeline_mng->getBoundaryRhsFE().reset();
pipeline_mng->getBoundaryLhsFE().reset();
auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
auto u_ptr = boost::make_shared<VectorDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"U", u_ptr}},
{},
{},
{})
);
pipeline_mng->getDomainRhsFE() = post_proc_fe;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR post_proc_fe->writeFile("out_result.h5m");
}
//! [Output results]
//! [Run program]
}
//! [Run program]
//! [Main]
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
MoFEM::Core::Initialize(&argc, &argv, param_file, help);
// Error handling
try {
// Register MoFEM discrete manager in PETSc
DMType dm_name = "DMMOFEM";
// Create MOAB instance
moab::Core mb_instance; // mesh database
moab::Interface &moab = mb_instance; // mesh database interface
// Create MoFEM instance
MoFEM::Core core(moab); // finite element database
MoFEM::Interface &m_field = core; // finite element interface
// Run the main analysis
Poisson2DiscontGalerkin poisson_problem(m_field);
CHKERR poisson_problem.runProgram();
}
// Finish work: cleaning memory, getting statistics, etc.
return 0;
}
//! [Main]
NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
Poisson2DiscontGalerkin::solveSystem
MoFEMErrorCode solveSystem()
[Set integration rules]
Definition: poisson_2d_dis_galerkin.cpp:236
Poisson2DiscontGalerkin::readMesh
MoFEMErrorCode readMesh()
[Read mesh]
Definition: poisson_2d_dis_galerkin.cpp:104
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
MoFEM::PipelineManager::getDomainRhsFE
boost::shared_ptr< FEMethod > & getDomainRhsFE()
Definition: PipelineManager.hpp:405
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
Poisson2DiscontGalerkin::runProgram
MoFEMErrorCode runProgram()
[Output results]
Definition: poisson_2d_dis_galerkin.cpp:520
H1
@ H1
continuous field
Definition: definitions.h:85
Poisson2DiscontGalerkin::checkResults
MoFEMErrorCode checkResults()
[Solve system]
Definition: poisson_2d_dis_galerkin.cpp:261
MoFEM::DataOperator
base operator to do operations at Gauss Pt. level
Definition: DataOperators.hpp:24
u_grad_exact
auto u_grad_exact
Definition: poisson_2d_dis_galerkin.cpp:49
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
DomainEle
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::DomainEle DomainEle
Definition: poisson_2d_dis_galerkin.cpp:16
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
phi
static double phi
Definition: poisson_2d_dis_galerkin.cpp:29
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
FaceSideEle
PipelineManager::ElementsAndOpsByDim< FE_DIM >::FaceSideEle FaceSideEle
Definition: level_set.cpp:41
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
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::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
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:523
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
Poisson2DiscontGalerkinOperators::OpL2LhsPenalty
Operator the left hand side matrix.
Definition: PoissonDiscontinousGalerkin.hpp:115
OpDomainGradGrad
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< BASE_DIM, FIELD_DIM, SPACE_DIM > OpDomainGradGrad
Definition: poisson_2d_dis_galerkin.cpp:36
DomainEleOp
DomainEle::UserDataOperator DomainEleOp
Definition: poisson_2d_dis_galerkin.cpp:17
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1576
FaceSideEle
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::FaceSideEle FaceSideEle
Definition: poisson_2d_dis_galerkin.cpp:23
Poisson2DiscontGalerkinOperators::OpL2BoundaryRhs
Operator to evaluate Dirichlet boundary conditions using DG.
Definition: PoissonDiscontinousGalerkin.hpp:274
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1294
Poisson2DiscontGalerkin::setupProblem
MoFEMErrorCode setupProblem()
[Read mesh]
Definition: poisson_2d_dis_galerkin.cpp:123
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
MoFEM::OpSetInvJacL2ForFace
Definition: UserDataOperators.hpp:2929
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
u_exact
auto u_exact
Definition: poisson_2d_dis_galerkin.cpp:42
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:667
MoFEM::Simple::getAddSkeletonFE
bool & getAddSkeletonFE()
Get the addSkeletonFE.
Definition: Simple.hpp:425
Poisson2DiscontGalerkin::simpleInterface
Simple * simpleInterface
Definition: poisson_2d_dis_galerkin.cpp:93
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:667
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
Poisson2DiscontGalerkin::domainField
std::string domainField
Definition: poisson_2d_dis_galerkin.cpp:96
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1067
Poisson2DiscontGalerkin::oRder
int oRder
Definition: poisson_2d_dis_galerkin.cpp:97
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
FaceSideOp
Poisson2DiscontGalerkin
Definition: poisson_2d_dis_galerkin.cpp:73
Poisson2DiscontGalerkinOperators
Definition: PoissonDiscontinousGalerkin.hpp:15
BoundaryEleOp
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
double
convert.type
type
Definition: convert.py:64
MoFEM::Problem::getNbDofsRow
DofIdx getNbDofsRow() const
Definition: ProblemsMultiIndices.hpp:376
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
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
Poisson2DiscontGalerkin::Poisson2DiscontGalerkin
Poisson2DiscontGalerkin(MoFEM::Interface &m_field)
Definition: poisson_2d_dis_galerkin.cpp:100
MoFEM::Simple::getAddBoundaryFE
bool & getAddBoundaryFE()
Get the addSkeletonFE.
Definition: Simple.hpp:436
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
SPACE_DIM
constexpr int SPACE_DIM
Definition: poisson_2d_dis_galerkin.cpp:12
Poisson2DiscontGalerkinOperators::OpCalculateSideData
Operator tp collect data from elements on the side of Edge/Face.
Definition: PoissonDiscontinousGalerkin.hpp:38
OpDomainSource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, FIELD_DIM > OpDomainSource
Definition: poisson_2d_dis_galerkin.cpp:38
nitsche
static double nitsche
Definition: poisson_2d_dis_galerkin.cpp:31
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
penalty
constexpr double penalty
Definition: shallow_wave.cpp:75
main
int main(int argc, char *argv[])
[Run program]
Definition: poisson_2d_dis_galerkin.cpp:537
BiLinearForm
BoundaryEle
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: poisson_2d_dis_galerkin.cpp:20
FIELD_DIM
constexpr int FIELD_DIM
Definition: poisson_2d_dis_galerkin.cpp:11
FTensor::Index< 'i', 2 >
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:473
PoissonDiscontinousGalerkin.hpp
help
static char help[]
Definition: poisson_2d_dis_galerkin.cpp:71
MoFEM::PipelineManager::getDomainLhsFE
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Definition: PipelineManager.hpp:401
DomainEleOp
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_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
OpDomainGradGrad
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
Definition: helmholtz.cpp:25
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
BoundaryEleOp
BoundaryEle::UserDataOperator BoundaryEleOp
Definition: poisson_2d_dis_galerkin.cpp:21
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1102
Poisson2DiscontGalerkin::boundaryCondition
MoFEMErrorCode boundaryCondition()
[Setup problem]
Definition: poisson_2d_dis_galerkin.cpp:165
Poisson2DiscontGalerkin::assembleSystem
MoFEMErrorCode assembleSystem()
[Boundary condition]
Definition: poisson_2d_dis_galerkin.cpp:173
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
Poisson2DiscontGalerkin::mField
MoFEM::Interface & mField
Definition: poisson_2d_dis_galerkin.cpp:92
MoFEM::DMMoFEMGetProblemPtr
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMoFEM.cpp:426
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition: PetscSmartObj.hpp:202
MoFEM::OpInvertMatrix
Definition: UserDataOperators.hpp:3254
MoFEM::OpCalculateHOJacForFace
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
Definition: HODataOperators.hpp:264
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
Poisson2DiscontGalerkin::outputResults
MoFEMErrorCode outputResults()
[Output results]
Definition: poisson_2d_dis_galerkin.cpp:477
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:54
is_test
PetscBool is_test
Definition: poisson_2d_dis_galerkin.cpp:40
EntData
EntitiesFieldData::EntData EntData
Definition: poisson_2d_dis_galerkin.cpp:14
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
source
auto source
Definition: poisson_2d_dis_galerkin.cpp:61
OpDomainSource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
Definition: child_and_parent.cpp:55
MoFEM::Simple::setUp
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:611
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
Poisson2DiscontGalerkin::setIntegrationRules
MoFEMErrorCode setIntegrationRules()
[Assemble system]
Definition: poisson_2d_dis_galerkin.cpp:215
BASE_DIM
constexpr int BASE_DIM
Definition: poisson_2d_dis_galerkin.cpp:10
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::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182