|
| v0.14.0
|
Go to the documentation of this file.
11 #ifndef __POISSON2DISGALERKIN_HPP__
12 #define __POISSON2DISGALERKIN_HPP__
22 std::array<VectorInt, 2>
24 std::array<VectorInt, 2>
43 std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE],
false);
45 for (
auto t = moab::CN::TypeDimensionMap[
SPACE_DIM].first;
46 t <= moab::CN::TypeDimensionMap[
SPACE_DIM].second; ++
t)
51 EntityType col_type,
EntData &row_data,
61 if ((CN::Dimension(row_type) ==
SPACE_DIM) &&
63 const auto nb_in_loop = getFEMethod()->nInTheLoop;
70 areaMap[nb_in_loop] = getMeasure();
71 senseMap[nb_in_loop] = getSkeletonSense();
92 &*base_mat.data().begin());
95 template <
typename T>
inline auto get_ntensor(T &base_mat,
int gg,
int bb) {
96 double *ptr = &base_mat(gg, bb);
101 double *ptr = &*base_mat.data().begin();
102 return getFTensor1FromPtr<2>(ptr);
105 template <
typename T>
107 double *ptr = &base_mat(gg, 2 * bb);
108 return getFTensor1FromPtr<2>(ptr);
132 const auto in_the_loop =
136 const std::array<std::string, 2> ele_type_name = {
"BOUNDARY",
"SKELETON"};
138 <<
"OpL2LhsPenalty inTheLoop " << ele_type_name[in_the_loop];
149 t_normal.normalize();
152 const size_t nb_integration_pts =
getGaussPts().size2();
156 const double beta =
static_cast<double>(
nitsche) / (in_the_loop + 1);
167 const auto sense_row =
senseMap[s0];
174 const auto sense_col =
senseMap[s1];
181 locMat.resize(nb_rows, nb_cols,
false);
190 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
195 auto t_mat =
locMat.data().begin();
199 for (; rr != nb_rows; ++rr) {
203 t_vn_plus(
i) = beta * (
phi * t_diff_row_base(
i) / p);
205 t_vn(
i) = t_row_base * t_normal(
i) * sense_row - t_vn_plus(
i);
209 auto t_diff_col_base =
213 for (
size_t cc = 0; cc != nb_cols; ++cc) {
217 t_un(
i) = -p * (t_col_base * t_normal(
i) * sense_col -
218 beta * t_diff_col_base(
i) / p);
221 *t_mat -= alpha * (t_vn(
i) * t_un(
i));
222 *t_mat -= alpha * (t_vn_plus(
i) * (beta * t_diff_col_base(
i)));
239 for (; rr < nb_row_base_functions; ++rr) {
252 &*
locMat.data().begin(), ADD_VALUES);
265 boost::shared_ptr<FaceSideEle>
285 const auto in_the_loop =
289 const double s = getMeasure() / (
areaMap[0]);
293 auto t_normal = getFTensor1Normal();
294 t_normal.normalize();
296 auto t_w = getFTensor0IntegrationWeight();
304 rhsF.resize(nb_rows,
false);
307 const size_t nb_integration_pts = getGaussPts().size2();
313 auto t_coords = getFTensor1CoordsAtGaussPts();
316 const double beta =
static_cast<double>(
nitsche) / (in_the_loop + 1);
318 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
319 const double alpha = getMeasure() * t_w;
321 const double source_val =
322 -p *
sourceFun(t_coords(0), t_coords(1), t_coords(2));
324 auto t_f =
rhsF.data().begin();
327 for (; rr != nb_rows; ++rr) {
330 t_vn_plus(
i) = beta * (
phi * t_diff_row_base(
i) / p);
332 t_vn(
i) = t_row_base * t_normal(
i) * sense_row - t_vn_plus(
i);
335 *t_f -= alpha * t_vn(
i) * (source_val * t_normal(
i));
342 for (; rr < nb_row_base_functions; ++rr) {
360 boost::shared_ptr<FaceSideEle>
368 #endif //__POISSON2DISGALERKIN_HPP__
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Data on single entity (This is passed as argument to DataOperator::doWork)
std::array< VectorInt, 2 > indicesRowSideMap
indices on rows for left hand-side
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
OpL2BoundaryRhs(boost::shared_ptr< FaceSideEle > side_fe_ptr, ScalarFun fun)
@ OPSPACE
operator do Work is execute on space data
OpL2LhsPenalty(boost::shared_ptr< FaceSideEle > side_fe_ptr)
Construct a new OpL2LhsPenalty.
VectorDouble rhsF
vector to store local operator right hand side
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
boost::shared_ptr< FaceSideEle > sideFEPtr
pointer to element to get data on edge/face sides
FTensor::Index< 'i', SPACE_DIM > i
std::array< MatrixDouble, 2 > rowBaseSideMap
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
auto get_diff_ntensor(T &base_mat)
Operator the left hand side matrix.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Operator to evaluate Dirichlet boundary conditions using DG.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
std::array< MatrixDouble, 2 > colDiffBaseSideMap
#define CHKERR
Inline error check.
boost::shared_ptr< FaceSideEle > sideFEPtr
pointer to element to get data on edge/face sides
auto getFTensor1Normal()
get normal as tensor
default operator for TRI element
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
const VectorInt & getIndices() const
Get global indices of dofs on entity.
auto get_ntensor(T &base_mat)
Operator tp collect data from elements on the side of Edge/Face.
std::array< MatrixDouble, 2 > rowDiffBaseSideMap
constexpr double t
plate stiffness
EntitiesFieldData::EntData EntData
constexpr auto field_name
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
#define MOFEM_LOG(channel, severity)
Log.
std::array< double, 2 > areaMap
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
std::array< MatrixDouble, 2 > colBaseSideMap
ScalarFun sourceFun
pointer to function to evaluate value of function on boundary
@ MOFEM_DATA_INCONSISTENCY
MatrixDouble locMat
local operator matrix
std::array< VectorInt, 2 > indicesColSideMap
indices on columns for left hand-side
UBlasVector< double > VectorDouble
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
auto fun
Function to approximate.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpCalculateSideData(std::string field_name, std::string col_field_name)
std::array< int, 2 > senseMap