v0.14.0
PoissonDiscontinousGalerkin.hpp

Example of implementation for discontinuous Galerkin.

/**
* \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] = getSkeletonSense();
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)
: BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE), sideFEPtr(side_fe_ptr) {}
MoFEMErrorCode doWork(int side, EntityType type,
// 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]);
auto t_w = getFTensor0IntegrationWeight();
// 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(),
&*locMat.data().begin(), ADD_VALUES);
// 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 Operator to evaluate Dirichlet boundary conditions using DG
*
*/
struct OpL2BoundaryRhs : public BoundaryEleOp {
public:
OpL2BoundaryRhs(boost::shared_ptr<FaceSideEle> side_fe_ptr, ScalarFun fun)
: BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE), sideFEPtr(side_fe_ptr),
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
// 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();
auto t_w = getFTensor0IntegrationWeight();
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(),
ADD_VALUES);
}
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 store local operator right hand side
};
}; // namespace Poisson2DiscontGalerkinOperators
#endif //__POISSON2DISGALERKIN_HPP__
NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
OpCalculateSideData
Operator tp collect data from elements on the side of Edge/Face.
Definition: plate.cpp:109
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
Poisson2DiscontGalerkinOperators::indicesRowSideMap
std::array< VectorInt, 2 > indicesRowSideMap
indices on rows for left hand-side
Definition: PoissonDiscontinousGalerkin.hpp:23
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1644
Poisson2DiscontGalerkinOperators::OpL2BoundaryRhs::OpL2BoundaryRhs
OpL2BoundaryRhs(boost::shared_ptr< FaceSideEle > side_fe_ptr, ScalarFun fun)
Definition: PoissonDiscontinousGalerkin.hpp:276
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
Poisson2DiscontGalerkinOperators::OpL2LhsPenalty::OpL2LhsPenalty
OpL2LhsPenalty(boost::shared_ptr< FaceSideEle > side_fe_ptr)
Construct a new OpL2LhsPenalty.
Definition: PoissonDiscontinousGalerkin.hpp:123
Poisson2DiscontGalerkinOperators::OpL2BoundaryRhs::rhsF
VectorDouble rhsF
vector to store local operator right hand side
Definition: PoissonDiscontinousGalerkin.hpp:363
Poisson2DiscontGalerkinOperators::LEFT_SIDE
@ LEFT_SIDE
Definition: PoissonDiscontinousGalerkin.hpp:20
phi
static double phi
Definition: poisson_2d_dis_galerkin.cpp:30
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
Poisson2DiscontGalerkinOperators::OpL2BoundaryRhs::sideFEPtr
boost::shared_ptr< FaceSideEle > sideFEPtr
pointer to element to get data on edge/face sides
Definition: PoissonDiscontinousGalerkin.hpp:361
Poisson2DiscontGalerkinOperators::i
FTensor::Index< 'i', SPACE_DIM > i
Definition: PoissonDiscontinousGalerkin.hpp:18
Poisson2DiscontGalerkinOperators::rowBaseSideMap
std::array< MatrixDouble, 2 > rowBaseSideMap
Definition: PoissonDiscontinousGalerkin.hpp:26
MoFEM::EntitiesFieldData::EntData::getDiffN
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
Definition: EntitiesFieldData.hpp:1329
Poisson2DiscontGalerkinOperators::get_diff_ntensor
auto get_diff_ntensor(T &base_mat)
Definition: PoissonDiscontinousGalerkin.hpp:100
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
FaceSideOp
FaceSideEle::UserDataOperator FaceSideOp
Definition: poisson_2d_dis_galerkin.cpp:25
Poisson2DiscontGalerkinOperators::colDiffBaseSideMap
std::array< MatrixDouble, 2 > colDiffBaseSideMap
Definition: PoissonDiscontinousGalerkin.hpp:29
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
FaceSideOp
Poisson2DiscontGalerkinOperators::OpL2LhsPenalty::sideFEPtr
boost::shared_ptr< FaceSideEle > sideFEPtr
pointer to element to get data on edge/face sides
Definition: PoissonDiscontinousGalerkin.hpp:266
Poisson2DiscontGalerkinOperators
Definition: PoissonDiscontinousGalerkin.hpp:15
BoundaryEleOp
convert.type
type
Definition: convert.py:64
Poisson2DiscontGalerkinOperators::OpCalculateSideData::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Definition: PoissonDiscontinousGalerkin.hpp:50
BoundaryEleOp
BoundaryEle::UserDataOperator BoundaryEleOp
Definition: child_and_parent.cpp:40
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1214
Poisson2DiscontGalerkinOperators::get_ntensor
auto get_ntensor(T &base_mat)
Definition: PoissonDiscontinousGalerkin.hpp:90
MoFEM::ScalarFun
boost::function< double(const double, const double, const double)> ScalarFun
Scalar function type.
Definition: FormsIntegrators.hpp:144
nitsche
static double nitsche
Definition: poisson_2d_dis_galerkin.cpp:32
Poisson2DiscontGalerkinOperators::rowDiffBaseSideMap
std::array< MatrixDouble, 2 > rowDiffBaseSideMap
Definition: PoissonDiscontinousGalerkin.hpp:28
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
penalty
constexpr double penalty
Definition: shallow_wave.cpp:75
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', SPACE_DIM >
Poisson2DiscontGalerkinOperators::OpL2BoundaryRhs::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Definition: PoissonDiscontinousGalerkin.hpp:280
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
Poisson2DiscontGalerkinOperators::areaMap
std::array< double, 2 > areaMap
Definition: PoissonDiscontinousGalerkin.hpp:30
FTensor::Tensor0
Definition: Tensor0.hpp:16
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: EntitiesFieldData.hpp:1318
Poisson2DiscontGalerkinOperators::colBaseSideMap
std::array< MatrixDouble, 2 > colBaseSideMap
Definition: PoissonDiscontinousGalerkin.hpp:27
Poisson2DiscontGalerkinOperators::OpL2BoundaryRhs::sourceFun
ScalarFun sourceFun
pointer to function to evaluate value of function on boundary
Definition: PoissonDiscontinousGalerkin.hpp:362
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
Poisson2DiscontGalerkinOperators::OpL2LhsPenalty::locMat
MatrixDouble locMat
local operator matrix
Definition: PoissonDiscontinousGalerkin.hpp:267
Poisson2DiscontGalerkinOperators::RIGHT_SIDE
@ RIGHT_SIDE
Definition: PoissonDiscontinousGalerkin.hpp:20
Poisson2DiscontGalerkinOperators::indicesColSideMap
std::array< VectorInt, 2 > indicesColSideMap
indices on columns for left hand-side
Definition: PoissonDiscontinousGalerkin.hpp:25
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
fun
auto fun
Function to approximate.
Definition: dg_projection.cpp:36
ElementSide
ElementSide
Definition: plate.cpp:92
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
Poisson2DiscontGalerkinOperators::OpL2LhsPenalty::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: PoissonDiscontinousGalerkin.hpp:126
Poisson2DiscontGalerkinOperators::OpCalculateSideData::OpCalculateSideData
OpCalculateSideData(std::string field_name, std::string col_field_name)
Definition: PoissonDiscontinousGalerkin.hpp:40
Poisson2DiscontGalerkinOperators::senseMap
std::array< int, 2 > senseMap
Definition: PoissonDiscontinousGalerkin.hpp:31