v0.14.0
Loading...
Searching...
No Matches
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 EntData = EntitiesFieldData::EntData;
using DomainEle = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::DomainEle;
using DomainEleOp = DomainEle::UserDataOperator;
using BoundaryEle =
PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
using BoundaryEleOp = BoundaryEle::UserDataOperator;
using FaceSideEle =
PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
using FaceSideOp = FaceSideEle::UserDataOperator;
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<
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
private:
// Declaration of other main functions called in runProgram()
// MoFEM interfaces
// Field name and approximation order
std::string domainField;
int oRder;
};
: domainField("U"), mField(m_field), oRder(4) {}
//! [Read mesh]
CHKERR simpleInterface->getOptions();
// 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.
simpleInterface->getAddSkeletonFE() = true;
simpleInterface->getAddBoundaryFE() = true;
CHKERR simpleInterface->loadFile();
}
//! [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);
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,
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;
CHKERR DMMoFEMGetProblemPtr(simpleInterface->getDM(), &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";
// 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]
static Index< 'o', 3 > o
static Index< 'p', 3 > p
std::string param_file
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
static char help[]
int main()
Definition: adol-c_atom.cpp:46
static const double eps
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
constexpr int SPACE_DIM
constexpr int FIELD_DIM
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
#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
static double penalty
ElementsAndOps< SPACE_DIM >::FaceSideEle FaceSideEle
static double phi
@ 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 DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMoFEM.cpp:412
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
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
constexpr int BASE_DIM
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)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
auto u_grad_exact
static double nitsche
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
BoundaryEle::UserDataOperator BoundaryEleOp
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)
@ 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.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
keeps basic data about problem
DofIdx getNbDofsRow() const
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Operator tp collect data from elements on the side of Edge/Face.
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 to evaluate Dirichlet boundary conditions using DG.