26#ifndef __DATAOPERATORS_HPP
27#define __DATAOPERATORS_HPP
29using namespace boost::numeric;
62 "doWork function is not implemented for this operator");
86 "doWork function is not implemented for this operator");
92 const bool error_if_no_base =
false);
96 std::array<bool, MBMAXTYPE>
132 template <
bool ErrorIfNoBase>
134 const std::array<bool, MBMAXTYPE> &do_entities);
150 :
tInvJac(&inv_jac(0, 0), &inv_jac(0, 1), &inv_jac(0, 2), &inv_jac(1, 0),
151 &inv_jac(1, 1), &inv_jac(1, 2), &inv_jac(2, 0), &inv_jac(2, 1),
173 :
tInvJac(&inv_jac(0, 0), &inv_jac(0, 1), &inv_jac(0, 2), &inv_jac(1, 0),
174 &inv_jac(1, 1), &inv_jac(1, 2), &inv_jac(2, 0), &inv_jac(2, 1),
207 tJac(&jac(0, 0), &jac(0, 1), &jac(0, 2), &jac(1, 0), &jac(1, 1),
208 &jac(1, 2), &jac(2, 0), &jac(2, 1), &jac(2, 2)) {}
238 :
tInvJac(&inv_jac(0, 0), &inv_jac(0, 1), &inv_jac(0, 2), &inv_jac(1, 0),
239 &inv_jac(1, 1), &inv_jac(1, 2), &inv_jac(2, 0), &inv_jac(2, 1),
273 template <
int R,
int D>
288 const int nb_base_functions = data.
getN().size2();
289 bool constant_diff =
false;
291 DIM * nb_base_functions) {
292 constant_diff =
true;
295 for (
unsigned int gg = 0; gg < data.
getN().size1(); gg++) {
296 double *data_ptr, *n_ptr, *diff_n_ptr;
297 n_ptr = &data.
getN()(gg, 0);
299 diff_n_ptr = &data.
getDiffN()(0, 0);
301 diff_n_ptr = &data.
getDiffN()(gg, 0);
304 for (
int rr = 0; rr < RANK; rr++, data_ptr++) {
306 cblas_ddot(nb_dofs / RANK, n_ptr, 1, data_ptr, RANK);
307 double *diff_n_ptr2 = diff_n_ptr;
308 for (
unsigned int dd = 0;
dd < DIM;
dd++, diff_n_ptr2++) {
310 cblas_ddot(nb_dofs / RANK, diff_n_ptr2, DIM, data_ptr, RANK);
330 if (nb_dofs % RANK != 0) {
332 "data inconsistency, type %d, side %d, nb_dofs %d, rank %d",
333 type, side, nb_dofs, RANK);
335 if (nb_dofs / RANK > data.
getN().size2()) {
336 std::cerr << side <<
" " <<
type <<
" "
338 std::cerr << data.
getN() << std::endl;
340 "data inconsistency nb_dofs >= data.N.size2(), i.e. %u >= %u",
341 nb_dofs, data.
getN().size2());
344 if (
type == MBVERTEX) {
363OpGetDataAndGradient<3, 3>::getValAtGaussPtsTensor<3>(
MatrixDouble &data);
371OpGetDataAndGradient<3, 3>::getGradAtGaussPtsTensor<3, 3>(
MatrixDouble &data);
447 const int normal_shift = 0)
454 const int normal_shift = 0)
#define MoFEMFunctionReturnHot(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 ...
@ MOFEM_DATA_INCONSISTENCY
static const char *const ApproximationBaseNames[]
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
implementation of Data Operators for Forces and Sources
base operator to do operations at Gauss Pt. level
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
virtual MoFEMErrorCode opLhs(EntitiesFieldData &row_data, EntitiesFieldData &col_data)
virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
bool & doPrisms
\deprectaed
bool sYmm
If true assume that matrix is symmetric structure.
boost::function< MoFEMErrorCode(DataOperator *op_ptr, int side, EntityType type, EntitiesFieldData::EntData &data)> DoWorkRhsHookFunType
boost::function< MoFEMErrorCode(DataOperator *op_ptr, int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)> DoWorkLhsHookFunType
void setSymm()
set if operator is executed taking in account symmetry
virtual MoFEMErrorCode opRhs(EntitiesFieldData &data, const bool error_if_no_base=false)
bool getSymm() const
Get if operator uses symmetry of DOFs or not.
virtual MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
bool & doVertices
\deprectaed If false skip vertices
DataOperator(const bool symm=true)
virtual ~DataOperator()=default
void unSetSymm()
unset if operator is executed for non symmetric problem
DoWorkLhsHookFunType doWorkLhsHook
DoWorkRhsHookFunType doWorkRhsHook
bool & doEdges
\deprectaed If false skip edges
bool & doQuads
\deprectaed
Data on single entity (This is passed as argument to DataOperator::doWork)
FieldApproximationBase & getBase()
Get approximation base.
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 VectorDouble & getFieldData() const
get dofs values
data structure for finite element entity
calculate normals at Gauss points of triangle element
MatrixDouble & tAngent2_at_GaussPtF3
MatrixDouble & tAngent1_at_GaussPtF3
MatrixDouble & nOrmals_at_GaussPtF3
MatrixDouble & cOords_at_GaussPtF3
MatrixDouble & nOrmals_at_GaussPtF4
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MatrixDouble & cOords_at_GaussPtF4
MoFEMErrorCode calculateNormals()
MatrixDouble & tAngent2_at_GaussPtF4
OpGetCoordsAndNormalsOnPrism(MatrixDouble &coords_at_gaussptf3, MatrixDouble &normals_at_gaussptf3, MatrixDouble &tangent1_at_gaussptf3, MatrixDouble &tangent2_at_gaussptf3, MatrixDouble &coords_at_gaussptf4, MatrixDouble &normals_at_gaussptf4, MatrixDouble &tangent1_at_gaussptf4, MatrixDouble &tangent2_at_gaussptf4)
MatrixDouble & tAngent1_at_GaussPtF4
Get field values and gradients at Gauss points.
FTensor::Tensor2< double *, R, D > getGradAtGaussPtsTensor(MatrixDouble &data)
MatrixDouble & dataGradAtGaussPts
MatrixDouble & dataAtGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode calculateValAndGrad(int side, EntityType type, EntitiesFieldData::EntData &data)
Calculate gradient and values at integration points.
OpGetDataAndGradient(MatrixDouble &data_at_gauss_pt, MatrixDouble &data_grad_at_gauss_pt)
FTensor::Tensor1< double *, R > getValAtGaussPtsTensor(MatrixDouble &data)
Calculate tangent vector on edge form HO geometry approximation.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpGetHOTangentOnEdge(MatrixDouble &tangent)
Transform local reference derivatives of shape function to global derivatives.
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpSetInvJacH1(MatrixDouble3by3 &inv_jac)
FTensor::Tensor2< double *, 3, 3 > tInvJac
brief Transform local reference derivatives of shape function to global derivatives
FTensor::Index< 'j', 3 > j
FTensor::Tensor2< double *, 3, 3 > tInvJac
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpSetInvJacHdivAndHcurl(MatrixDouble3by3 &inv_jac)
FTensor::Index< 'i', 3 > i
FTensor::Index< 'k', 3 > k
MatrixDouble diffHdivInvJac