v0.13.1
SCL-11: Discountinous Galrekin for Poisson problem

Note
Prerequisites of this tutorial include SCL-1: Poisson's equation (homogeneous BC)

Note
Intended learning outcome:

# Problem formulation

$-(v_{,j}, \sigma_{j}(\mathbf{u}))_V+(v, \overline{t})_{S^\sigma} = 0$

$-(v_{j}, \sigma_{j}(\mathbf{u}))_\Omega -\sum_{i} (\left\{[v n_j] - \theta\gamma\sigma_{j}(\mathbf{v})\right\} + \theta\gamma\sigma_{j}(\mathbf{v}), \hat{\sigma}_{j}(\mathbf{u}))_{F_i}+(v_i, \overline{t}_i)_{S^\sigma}= 0$

$-(v_{j}, \sigma_{j}(\mathbf{u}))_\Omega+(v_i, \overline{t})_{S^\sigma} -\sum_{i} (\left\{[v_i n_j] - \theta\gamma\sigma_{ij}(\mathbf{v})\right\}, \hat{\sigma}_{j}(\mathbf{u}))_{F_i}-\sum_{i} (\theta\gamma\sigma_{j}(\mathbf{v}),\sigma_{j}(\mathbf{u}))_{F_i} = 0$

$\hat{\sigma} = \gamma^{-1} \left\{ [u n_j] - \gamma \sigma_{j} (\mathbf{u}) \right\}$

$-(v_{j}, \sigma_{j}(\mathbf{u}))_\Omega+(v, \overline{t})_{S^\sigma} -\sum_{i} (\left\{[v n_j] - \theta\gamma\sigma_{j}(\mathbf{v})\right\}, \gamma^{-1} \left\{ [u n_j] - \gamma \sigma_{j}(\mathbf{u})\right\})_{F_i} -\sum_{i} (\theta\gamma\sigma_{j}(\mathbf{v}), \gamma^{-1} \left\{\gamma \sigma_{j}(\mathbf{u})\right\})_{F_i} = 0$

# Source code main code

/**
* \file poisson_2d_dis_galerkin.cpp
* \example poisson_2d_dis_galerkin.cpp
*
* Example of implementation for discontinuous Galerkin.
*/
#include <BasicFiniteElements.hpp>
template <int DIM> struct ElementsAndOps {};
template <> struct ElementsAndOps<2> {
using FaceSideEle = FaceElementForcesAndSourcesCoreOnSide;
};
constexpr int BASE_DIM = 1;
constexpr int FIELD_DIM = 1;
constexpr int SPACE_DIM = 2;
using PostProcEle = PostProcBrokenMeshInMoab<DomainEle>;
static double penalty = 1e6;
static double phi =
-1; // 1 - symmetric Nitsche, 0 - nonsymmetric, -1 antisymmetrica
static double nitsche = 1;
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) {}
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.
}
//! [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));
};
[](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);
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>();
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_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));
error[L2] += alpha * pow(diff, 2);
++t_w;
++t_val;
++t_coords;
}
error[H1] = error[L2] + error[SEMI_NORM];
CHKERR VecSetValues(l2_vec, LAST_NORM, error_indices.data(), error.data(),
}
};
auto side_fe_ptr = boost::make_shared<FaceSideEle>(mField);
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(),
};
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;
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");
}
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
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:304
static char help[]
MoFEM::FaceElementForcesAndSourcesCore FaceEle
static const double eps
EntitiesFieldData::EntData EntData
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
constexpr int SPACE_DIM
constexpr int FIELD_DIM
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
BoundaryEle::UserDataOperator BoundaryEleOp
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
int main(int argc, char *argv[])
#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
static double phi
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
constexpr int BASE_DIM
FTensor::Index< 'i', SPACE_DIM > i
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
static double nitsche
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: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.
Definition: PostProc.hpp:166
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
DofIdx getNbDofsRow() const
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
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()
MoFEMErrorCode setIntegrationRules()
[Assemble system]
MoFEMErrorCode assembleSystem()
[Boundary condition]
Poisson2DiscontGalerkin(MoFEM::Interface &m_field)
MoFEMErrorCode solveSystem()
[Set integration rules]
MoFEMErrorCode outputResults()
[Output results]
MoFEMErrorCode runProgram()
[Output results]
Opator tp evaluate Dirichlet boundary conditions using DG.

# Source code main code

/**
* \file PoissonDiscontinousGalerkin.hpp
* \example PoissonDiscontinousGalerkin.hpp
*
* Example of implementation for discontinuous Galerkin.
*/
// Define name if it has not been defined yet
#ifndef __POISSON2DISGALERKIN_HPP__
#define __POISSON2DISGALERKIN_HPP__
// Namespace that contains necessary UDOs, will be included in the main program
// Declare FTensor index for 2D problem
// data for skeleton computation
std::array<VectorInt, 2>
indicesRowSideMap; ///< indices on rows for left hand-side
std::array<VectorInt, 2>
indicesColSideMap; ///< indices on columns for left hand-side
std::array<MatrixDouble, 2> rowBaseSideMap; // base functions on rows
std::array<MatrixDouble, 2> colBaseSideMap; // base function on columns
std::array<MatrixDouble, 2> rowDiffBaseSideMap; // derivative of base functions
std::array<MatrixDouble, 2> colDiffBaseSideMap; // derivative of base functions
std::array<double, 2> areaMap; // area/volume of elements on the side
std::array<int, 2> senseMap; // orientaton of local element edge/face in respect
// to global orientation of edge/face
/**
* @brief Operator tp collect data from elements on the side of Edge/Face
*
*/
OpCalculateSideData(std::string field_name, std::string col_field_name)
: FaceSideOp(field_name, col_field_name, FaceSideOp::OPROWCOL) {
std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
for (auto t = moab::CN::TypeDimensionMap[SPACE_DIM].first;
t <= moab::CN::TypeDimensionMap[SPACE_DIM].second; ++t)
doEntities[t] = true;
}
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type, EntData &row_data,
EntData &col_data) {
// Note: THat for L2 base data rows, and columns are the same, so operator
// above can be simpler operator for the right hand side, and data can be
// stored only for rows, since for columns data are the same. However for
// complex multi-physics problems that not necessary would be a case. For
// generality, we keep generic case.
if ((CN::Dimension(row_type) == SPACE_DIM) &&
(CN::Dimension(col_type) == SPACE_DIM)) {
const auto nb_in_loop = getFEMethod()->nInTheLoop;
indicesRowSideMap[nb_in_loop] = row_data.getIndices();
indicesColSideMap[nb_in_loop] = col_data.getIndices();
rowBaseSideMap[nb_in_loop] = row_data.getN();
colBaseSideMap[nb_in_loop] = col_data.getN();
rowDiffBaseSideMap[nb_in_loop] = row_data.getDiffN();
colDiffBaseSideMap[nb_in_loop] = col_data.getDiffN();
areaMap[nb_in_loop] = getMeasure();
senseMap[nb_in_loop] = getEdgeSense();
if (!nb_in_loop) {
indicesRowSideMap[1].clear();
indicesColSideMap[1].clear();
rowBaseSideMap[1].clear();
colBaseSideMap[1].clear();
rowDiffBaseSideMap[1].clear();
colDiffBaseSideMap[1].clear();
areaMap[1] = 0;
senseMap[1] = 0;
}
} else {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Should not happen");
}
}
};
template <typename T> inline auto get_ntensor(T &base_mat) {
&*base_mat.data().begin());
};
template <typename T> inline auto get_ntensor(T &base_mat, int gg, int bb) {
double *ptr = &base_mat(gg, bb);
};
template <typename T> inline auto get_diff_ntensor(T &base_mat) {
double *ptr = &*base_mat.data().begin();
return getFTensor1FromPtr<2>(ptr);
};
template <typename T>
inline auto get_diff_ntensor(T &base_mat, int gg, int bb) {
double *ptr = &base_mat(gg, 2 * bb);
return getFTensor1FromPtr<2>(ptr);
};
/**
* @brief Operator the left hand side matrix
*
*/
struct OpL2LhsPenalty : public BoundaryEleOp {
public:
/**
* @brief Construct a new OpL2LhsPenalty
*
* @param side_fe_ptr pointer to FE to evaluate side elements
*/
OpL2LhsPenalty(boost::shared_ptr<FaceSideEle> side_fe_ptr)
// Collect data from side domian elements
CHKERR loopSideFaces("dFE", *sideFEPtr);
const auto in_the_loop =
sideFEPtr->nInTheLoop; // return number of elements on the side
#ifndef NDEBUG
const std::array<std::string, 2> ele_type_name = {"BOUNDARY", "SKELETON"};
MOFEM_LOG("SELF", Sev::noisy)
<< "OpL2LhsPenalty inTheLoop " << ele_type_name[in_the_loop];
MOFEM_LOG("SELF", Sev::noisy)
<< "OpL2LhsPenalty sense " << senseMap[0] << " " << senseMap[1];
#endif
// calculate penalty
const double s = getMeasure() / (areaMap[0] + areaMap[1]);
const double p = penalty * s;
// get normal of the face or edge
auto t_normal = getFTensor1Normal();
t_normal.normalize();
// get number of integration points on face
const size_t nb_integration_pts = getGaussPts().size2();
// beta paramter is zero, when penalty method is used, also, takes value 1,
// when boundary edge/face is evaluated, and 2 when is skeleton edge/face.
const double beta = static_cast<double>(nitsche) / (in_the_loop + 1);
// iterate the sides rows
for (auto s0 : {LEFT_SIDE, RIGHT_SIDE}) {
// gent number of DOFs on the right side.
const auto nb_rows = indicesRowSideMap[s0].size();
if (nb_rows) {
// get orientation of the local element edge
const auto sense_row = senseMap[s0];
// iterate the side cols
const auto nb_row_base_functions = rowBaseSideMap[s0].size2();
for (auto s1 : {LEFT_SIDE, RIGHT_SIDE}) {
// get orientation of the edge
const auto sense_col = senseMap[s1];
// number of dofs, for homogenous approximation this value is
// constant.
const auto nb_cols = indicesColSideMap[s1].size();
// resize local element matrix
locMat.resize(nb_rows, nb_cols, false);
locMat.clear();
// get base functions, and integration weights
auto t_row_base = get_ntensor(rowBaseSideMap[s0]);
auto t_diff_row_base = get_diff_ntensor(rowDiffBaseSideMap[s0]);
// iterate integration points on face/edge
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
// t_w is integration weight, and measure is element area, or
// volume, depending if problem is in 2d/3d.
const double alpha = getMeasure() * t_w;
auto t_mat = locMat.data().begin();
// iterate rows
size_t rr = 0;
for (; rr != nb_rows; ++rr) {
// calculate tetting function
t_vn_plus(i) = beta * (phi * t_diff_row_base(i) / p);
t_vn(i) = t_row_base * t_normal(i) * sense_row - t_vn_plus(i);
// get base functions on columns
auto t_col_base = get_ntensor(colBaseSideMap[s1], gg, 0);
auto t_diff_col_base =
// iterate columns
for (size_t cc = 0; cc != nb_cols; ++cc) {
// calculate variance of tested function
t_un(i) = -p * (t_col_base * t_normal(i) * sense_col -
beta * t_diff_col_base(i) / p);
// assemble matrix
*t_mat -= alpha * (t_vn(i) * t_un(i));
*t_mat -= alpha * (t_vn_plus(i) * (beta * t_diff_col_base(i)));
// move to next column base and element of matrix
++t_col_base;
++t_diff_col_base;
++t_mat;
}
// move to next row base
++t_row_base;
++t_diff_row_base;
}
// this is obsolete for this particular example, we keep it for
// generality. in case of multi-physics problems different fields can
// chare hierarchical base but use different approx. order, so is
// possible to have more base functions than DOFs on element.
for (; rr < nb_row_base_functions; ++rr) {
++t_row_base;
++t_diff_row_base;
}
++t_w;
}
// assemble system
&*indicesRowSideMap[s0].begin(),
indicesColSideMap[s1].size(),
&*indicesColSideMap[s1].begin(),
// stop of boundary element
if (!in_the_loop)
}
}
}
}
private:
boost::shared_ptr<FaceSideEle>
sideFEPtr; ///< pointer to element to get data on edge/face sides
MatrixDouble locMat; ///< local operator matrix
};
/**
* @brief Opator tp evaluate Dirichlet boundary conditions using DG
*
*/
struct OpL2BoundaryRhs : public BoundaryEleOp {
public:
OpL2BoundaryRhs(boost::shared_ptr<FaceSideEle> side_fe_ptr, ScalarFun fun)
// get normal of the face or edge
CHKERR loopSideFaces("dFE", *sideFEPtr);
const auto in_the_loop =
sideFEPtr->nInTheLoop; // return number of elements on the side
// calculate penalty
const double s = getMeasure() / (areaMap[0]);
const double p = penalty * s;
// get normal of the face or edge
auto t_normal = getFTensor1Normal();
t_normal.normalize();
const size_t nb_rows = indicesRowSideMap[0].size();
if (!nb_rows)
// resize local operator vector
rhsF.resize(nb_rows, false);
rhsF.clear();
const size_t nb_integration_pts = getGaussPts().size2();
const size_t nb_row_base_functions = rowBaseSideMap[0].size2();
// shape funcs
auto t_row_base = get_ntensor(rowBaseSideMap[0]);
auto t_diff_row_base = get_diff_ntensor(rowDiffBaseSideMap[0]);
auto t_coords = getFTensor1CoordsAtGaussPts();
const auto sense_row = senseMap[0];
const double beta = static_cast<double>(nitsche) / (in_the_loop + 1);
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = getMeasure() * t_w;
const double source_val =
-p * sourceFun(t_coords(0), t_coords(1), t_coords(2));
auto t_f = rhsF.data().begin();
size_t rr = 0;
for (; rr != nb_rows; ++rr) {
t_vn_plus(i) = beta * (phi * t_diff_row_base(i) / p);
t_vn(i) = t_row_base * t_normal(i) * sense_row - t_vn_plus(i);
// assemble operator local vactor
*t_f -= alpha * t_vn(i) * (source_val * t_normal(i));
++t_row_base;
++t_diff_row_base;
++t_f;
}
for (; rr < nb_row_base_functions; ++rr) {
++t_row_base;
++t_diff_row_base;
}
++t_w;
++t_coords;
}
// assemble local operator vector to global vector
&*indicesRowSideMap[0].begin(), &*rhsF.data().begin(),
}
private:
boost::shared_ptr<FaceSideEle>
sideFEPtr; ///< pointer to element to get data on edge/face sides
ScalarFun sourceFun; ///< pointer to function to evaluate value of function on boundary
VectorDouble rhsF; ///< vector to strore local operator right hand side
};
}; // namespace Poisson2DiscontGalerkinOperators
#endif //__POISSON2DISGALERKIN_HPP__
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
boost::function< double(const double, const double, const double)> ScalarFun
Scalar function type.
auto fun
Function to approximate.
const double T
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
FTensor::Tensor1< FTensor::PackPtr< double *, 2 >, 2 > getFTensor1FromPtr< 2 >(double *ptr)
Definition: Templates.hpp:726
std::array< MatrixDouble, 2 > rowBaseSideMap
std::array< VectorInt, 2 > indicesColSideMap
indices on columns for left hand-side
FTensor::Index< 'i', SPACE_DIM > i
std::array< MatrixDouble, 2 > colDiffBaseSideMap
std::array< MatrixDouble, 2 > colBaseSideMap
std::array< VectorInt, 2 > indicesRowSideMap
indices on rows for left hand-side
std::array< MatrixDouble, 2 > rowDiffBaseSideMap
constexpr double t
plate stiffness
Definition: plate.cpp:59
constexpr auto field_name
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of dofs on entity.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
auto getFTensor0IntegrationWeight()
Get integration weights.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
OpCalculateSideData(std::string field_name, std::string col_field_name)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
VectorDouble rhsF
vector to strore local operator right hand side
OpL2BoundaryRhs(boost::shared_ptr< FaceSideEle > side_fe_ptr, ScalarFun fun)
ScalarFun sourceFun
pointer to function to evaluate value of function on boundary
boost::shared_ptr< FaceSideEle > sideFEPtr
pointer to element to get data on edge/face sides
OpL2LhsPenalty(boost::shared_ptr< FaceSideEle > side_fe_ptr)
Construct a new OpL2LhsPenalty.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
boost::shared_ptr< FaceSideEle > sideFEPtr
pointer to element to get data on edge/face sides