14#ifndef __DATAOPERATORS_HPP
15#define __DATAOPERATORS_HPP
17using namespace boost::numeric;
50 "doWork function is not implemented for this operator");
74 "doWork function is not implemented for this operator");
80 const bool error_if_no_base =
false);
84 std::array<bool, MBMAXTYPE>
120 template <
bool ErrorIfNoBase>
122 const std::array<bool, MBMAXTYPE> &do_entities);
138 :
tInvJac(&inv_jac(0, 0), &inv_jac(0, 1), &inv_jac(0, 2), &inv_jac(1, 0),
139 &inv_jac(1, 1), &inv_jac(1, 2), &inv_jac(2, 0), &inv_jac(2, 1),
161 :
tInvJac(&inv_jac(0, 0), &inv_jac(0, 1), &inv_jac(0, 2), &inv_jac(1, 0),
162 &inv_jac(1, 1), &inv_jac(1, 2), &inv_jac(2, 0), &inv_jac(2, 1),
195 tJac(&jac(0, 0), &jac(0, 1), &jac(0, 2), &jac(1, 0), &jac(1, 1),
196 &jac(1, 2), &jac(2, 0), &jac(2, 1), &jac(2, 2)) {}
226 :
tInvJac(&inv_jac(0, 0), &inv_jac(0, 1), &inv_jac(0, 2), &inv_jac(1, 0),
227 &inv_jac(1, 1), &inv_jac(1, 2), &inv_jac(2, 0), &inv_jac(2, 1),
261 template <
int R,
int D>
276 const int nb_base_functions = data.
getN().size2();
277 bool constant_diff =
false;
278 if (type == MBVERTEX && data.
getDiffN().size1() * data.
getDiffN().size2() ==
279 DIM * nb_base_functions) {
280 constant_diff =
true;
283 for (
unsigned int gg = 0; gg < data.
getN().size1(); gg++) {
284 double *data_ptr, *n_ptr, *diff_n_ptr;
285 n_ptr = &data.
getN()(gg, 0);
287 diff_n_ptr = &data.
getDiffN()(0, 0);
289 diff_n_ptr = &data.
getDiffN()(gg, 0);
292 for (
int rr = 0; rr < RANK; rr++, data_ptr++) {
294 cblas_ddot(nb_dofs / RANK, n_ptr, 1, data_ptr, RANK);
295 double *diff_n_ptr2 = diff_n_ptr;
296 for (
unsigned int dd = 0; dd < DIM; dd++, diff_n_ptr2++) {
298 cblas_ddot(nb_dofs / RANK, diff_n_ptr2, DIM, data_ptr, RANK);
318 if (nb_dofs % RANK != 0) {
320 "data inconsistency, type %d, side %d, nb_dofs %d, rank %d",
321 type, side, nb_dofs, RANK);
323 if (nb_dofs / RANK > data.
getN().size2()) {
324 std::cerr << side <<
" " << type <<
" "
326 std::cerr << data.
getN() << std::endl;
328 "data inconsistency nb_dofs >= data.N.size2(), i.e. %u >= %u",
329 nb_dofs, data.
getN().size2());
332 if (type == MBVERTEX) {
351OpGetDataAndGradient<3, 3>::getValAtGaussPtsTensor<3>(
MatrixDouble &data);
359OpGetDataAndGradient<3, 3>::getGradAtGaussPtsTensor<3, 3>(
MatrixDouble &data);
435 const int normal_shift = 0)
442 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.
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
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