v0.13.1
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.
*/
/* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
#include <BasicFiniteElements.hpp>
template <int DIM> struct ElementsAndOps {};
template <> struct ElementsAndOps<2> {
};
constexpr int BASE_DIM = 1;
constexpr int FIELD_DIM = 1;
constexpr int SPACE_DIM = 2;
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 diminesion (edge) elements, despite that fact that
// there is no field spanning on such elements. We need them for DF 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 = smartCreateDMVector(dm);
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 = createSmartVectorMPI(
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);
post_proc_fe->generateReferenceElementMesh();
post_proc_fe->addFieldValuesPostProc(domainField);
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
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:314
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:385
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:75
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ H1
continuous field
Definition: definitions.h:98
@ NOSPACE
Definition: definitions.h:96
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:53
#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
constexpr double alpha
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:482
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:385
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:977
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:311
FTensor::Index< 'i', SPACE_DIM > i
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, 2 > OpDomainGradGrad
Definition: helmholtz.cpp:37
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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:1096
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
const double D
diffusivity
int main(int argc, char *argv[])
[Run program]
EntitiesFieldData::EntData EntData
auto u_grad_exact
ElementsAndOps< SPACE_DIM >::DomainEleOp DomainEleOp
static char help[]
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::BoundaryEleOp BoundaryEleOp
static double nitsche
constexpr int FIELD_DIM
constexpr int BASE_DIM
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< BASE_DIM, FIELD_DIM, SPACE_DIM > OpDomainGradGrad
ElementsAndOps< SPACE_DIM >::FaceSideEle FaceSideEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, FIELD_DIM > OpDomainSource
static double phi
PetscBool is_test
constexpr double penalty
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: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
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.
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:35
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Operator tp collect data from elements on the side of Edge/Face.
Definition: plate.cpp:127
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]
Opator tp evaluate Dirichlet boundary conditions using DG.
Postprocess on face.