7#ifndef __POISSONOPERATORS_HPP__ 
    8#define __POISSONOPERATORS_HPP__ 
   62    if (row_side == col_side && row_type == col_type) {
 
 
   99    double vol = getVolume();
 
  101    auto t_w = getFTensor0IntegrationWeight();
 
  107      const double alpha = t_w * vol;
 
  113      for (
int rr = 0; rr != 
nbRows; rr++) {
 
  117        for (
int cc = 0; cc != 
nbCols; cc++) {
 
  119          a += alpha * (t_row_grad(
i) * t_col_grad(
i));
 
 
  141    const int *row_indices = &*row_data.
getIndices().data().begin();
 
  143    const int *col_indices = &*col_data.
getIndices().data().begin();
 
  144    Mat 
B = getFEMethod()->ksp_B != PETSC_NULLPTR ? getFEMethod()->ksp_B
 
  145                                               : getFEMethod()->snes_B;
 
  148                        &*
locMat.data().begin(), ADD_VALUES);
 
  155                          &*
locMat.data().begin(), ADD_VALUES);
 
 
 
  220    : 
public OpBaseRhs<VolumeElementForcesAndSourcesCore::UserDataOperator> {
 
  222  typedef boost::function<
double(
const double, 
const double, 
const double)>
 
  249    double vol = getVolume();
 
  251    auto t_w = getFTensor0IntegrationWeight();
 
  255    auto t_coords = getFTensor1CoordsAtGaussPts();
 
  260          vol * t_w * 
fSource(t_coords(
NX), t_coords(
NY), t_coords(
NZ));
 
  265      for (
int rr = 0; rr != 
nbRows; rr++) {
 
 
  285    const int *indices = &*data.
getIndices().data().begin();
 
  287    const double *vals = &*
locVec.data().begin();
 
  288    Vec f = getFEMethod()->ksp_f != PETSC_NULLPTR ? getFEMethod()->ksp_f
 
  289                                               : getFEMethod()->snes_f;
 
 
 
  307  OpC(
const bool assemble_transpose)
 
 
  356    const double area = getArea();
 
  358    auto t_w = getFTensor0IntegrationWeight();
 
  363      const double alpha = area * t_w;
 
  368      for (
int rr = 0; rr != 
nbRows; rr++) {
 
  372        for (
int cc = 0; cc != 
nbCols; cc++) {
 
  374          c += alpha * t_row * t_col;
 
 
  392    const int *row_indices = &*row_data.
getIndices().data().begin();
 
  394    const int *col_indices = &*col_data.
getIndices().data().begin();
 
  395    Mat 
B = getFEMethod()->ksp_B != PETSC_NULLPTR ? getFEMethod()->ksp_B
 
  396                                               : getFEMethod()->snes_B;
 
  399                        &*
locMat.data().begin(), ADD_VALUES);
 
  405                          &*
locMat.data().begin(), ADD_VALUES);
 
 
 
  420    : 
public OpBaseRhs<FaceElementForcesAndSourcesCore::UserDataOperator> {
 
  422  typedef boost::function<
double(
const double, 
const double, 
const double)>
 
  449    const double area = getArea() * 
bEta;
 
  451    auto t_w = getFTensor0IntegrationWeight();
 
  455    auto t_coords = getFTensor1CoordsAtGaussPts();
 
  461          area * t_w * 
fValue(t_coords(
NX), t_coords(
NY), t_coords(
NZ));
 
  466      for (
int rr = 0; rr != 
nbRows; rr++) {
 
 
  482    const int *indices = &*data.
getIndices().data().begin();
 
  483    const double *vals = &*
locVec.data().begin();
 
  484    Vec f = getFEMethod()->ksp_f != PETSC_NULLPTR ? getFEMethod()->ksp_f
 
  485                                               : getFEMethod()->snes_f;
 
 
 
  495    : 
public OpBaseRhs<VolumeElementForcesAndSourcesCore::UserDataOperator> {
 
  497  typedef boost::function<
double(
const double, 
const double, 
const double)>
 
  499  typedef boost::function<FTensor::Tensor1<double, 3>(
 
  500      const double, 
const double, 
const double)>
 
  504          boost::shared_ptr<VectorDouble> &field_vals,
 
  505          boost::shared_ptr<MatrixDouble> &grad_vals, Vec global_error)
 
 
  543    const double vol = getVolume();
 
  545    auto t_w = getFTensor0IntegrationWeight();
 
  549    auto t_grad = getFTensor1FromMat<3>(*
gradVals);
 
  551    auto t_coords = getFTensor1CoordsAtGaussPts();
 
  556      double alpha = vol * t_w;
 
  558      double exact_u = 
uValue(t_coords(
NX), t_coords(
NY), t_coords(
NZ));
 
  560      t_exact_grad = 
gValue(t_coords(
NX), t_coords(
NY), t_coords(
NZ));
 
  562      t_error_grad(
i) = t_grad(
i) - t_exact_grad(
i);
 
  564      double error = pow(t_u - exact_u, 2) + t_error_grad(
i) * t_error_grad(
i);
 
 
 
  590  OpKt(boost::function<
double(
const double)> 
a,
 
  591       boost::function<
double(
const double)> diff_a,
 
  592       boost::shared_ptr<VectorDouble> &field_vals,
 
  593       boost::shared_ptr<MatrixDouble> &grad_vals)
 
 
  612    double vol = getVolume();
 
  614    auto t_w = getFTensor0IntegrationWeight();
 
  618    auto t_grad = getFTensor1FromMat<3>(*
gradVals);
 
  624      const double alpha = t_w * vol;
 
  625      const double beta = alpha * 
A(t_u);
 
  627      t_gamma(
i) = (alpha * 
diffA(t_u)) * t_grad(
i);
 
  632      for (
int rr = 0; rr != 
nbRows; rr++) {
 
  638        for (
int cc = 0; cc != 
nbCols; cc++) {
 
  640          a += (t_row_grad(
i) * beta) * t_col_grad(
i) +
 
  641               t_row_grad(
i) * (t_gamma(
i) * t_col);
 
 
 
  665                boost::shared_ptr<VectorDouble> &field_vals,
 
  666                boost::shared_ptr<MatrixDouble> &grad_vals)
 
 
  684    double vol = getVolume();
 
  686    auto t_w = getFTensor0IntegrationWeight();
 
  690    auto t_grad = getFTensor1FromMat<3>(*
gradVals);
 
  696    auto t_coords = getFTensor1CoordsAtGaussPts();
 
  700      const double alpha = vol * t_w;
 
  701      const double source_term =
 
  704      grad_term(
i) = (alpha * 
A(t_u)) * t_grad(
i);
 
  709      for (
int rr = 0; rr != 
nbRows; rr++) {
 
  711        t_a += t_v_grad(
i) * grad_term(
i) + t_v * source_term;
 
 
 
  731  OpRes_g(
FVal f_value, boost::shared_ptr<VectorDouble> &field_vals)
 
 
  747    const double area = getArea() * 
bEta;
 
  749    auto t_w = getFTensor0IntegrationWeight();
 
  755    auto t_coords = getFTensor1CoordsAtGaussPts();
 
  760      double alpha = area * t_w;
 
  764      for (
int rr = 0; rr != 
nbRows; rr++) {
 
  766               (t_u - 
fValue(t_coords(
NX), t_coords(
NY), t_coords(
NZ)));
 
 
 
  796    const double area = getArea() * 
bEta;
 
  798    auto t_w = getFTensor0IntegrationWeight();
 
  807      double alpha = area * t_w;
 
  811      for (
int rr = 0; rr != 
nbRows; rr++) {
 
  812        t_a += alpha * t_u * t_lambda;
 
 
 
  835  int operator()(
int, 
int, 
int p)
 const { 
return 2 * (p - 1); }
 
 
  851    return 2 * p_data + 1;
 
 
 
  869      boost::function<
double(
const double, 
const double, 
const double)> f_u,
 
  870      boost::function<
double(
const double, 
const double, 
const double)>
 
  872      boost::shared_ptr<ForcesAndSourcesCore> &domain_lhs_fe,
 
  873      boost::shared_ptr<ForcesAndSourcesCore> &boundary_lhs_fe,
 
  874      boost::shared_ptr<ForcesAndSourcesCore> &domain_rhs_fe,
 
  875      boost::shared_ptr<ForcesAndSourcesCore> &boundary_rhs_fe,
 
  876      bool trans = 
true)
 const {
 
  880    domain_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
 
  882    boundary_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
 
  884    domain_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
 
  886    boundary_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
 
  890    domain_lhs_fe->getRuleHook = 
VolRule();
 
  891    domain_rhs_fe->getRuleHook = 
VolRule();
 
  892    boundary_lhs_fe->getRuleHook = 
FaceRule();
 
  893    boundary_rhs_fe->getRuleHook = 
FaceRule();
 
  897    domain_lhs_fe->getOpPtrVector().push_back(
new OpK());
 
  899    domain_rhs_fe->getOpPtrVector().push_back(
new OpF(f_source));
 
  901    boundary_lhs_fe->getOpPtrVector().push_back(
new OpC(trans));
 
  903    boundary_rhs_fe->getOpPtrVector().push_back(
new Op_g(f_u));
 
 
  912      boost::function<
double(
const double, 
const double, 
const double)> f_u,
 
  917      boost::shared_ptr<ForcesAndSourcesCore> &domain_error)
 const {
 
  920    domain_error = boost::shared_ptr<ForcesAndSourcesCore>(
 
  922    domain_error->getRuleHook = 
VolRule();
 
  926    boost::shared_ptr<VectorDouble> values_at_integration_ptr =
 
  927        boost::make_shared<VectorDouble>();
 
  929    boost::shared_ptr<MatrixDouble> grad_at_integration_ptr =
 
  930        boost::make_shared<MatrixDouble>();
 
  932    domain_error->getOpPtrVector().push_back(
 
  935    domain_error->getOpPtrVector().push_back(
 
  938    domain_error->getOpPtrVector().push_back(
 
  939        new OpError(f_u, g_u, values_at_integration_ptr,
 
  940                    grad_at_integration_ptr, global_error));
 
 
  948      boost::shared_ptr<PostProcFE> &post_proc_volume)
 const {
 
  965    auto det_ptr = boost::make_shared<VectorDouble>();
 
  966    auto jac_ptr = boost::make_shared<MatrixDouble>();
 
  967    auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
 
  968    post_proc_volume->getOpPtrVector().push_back(
 
  970    post_proc_volume->getOpPtrVector().push_back(
 
  972    post_proc_volume->getOpPtrVector().push_back(
 
  975    auto u_ptr = boost::make_shared<VectorDouble>();
 
  976    auto grad_ptr = boost::make_shared<MatrixDouble>();
 
  977    auto e_ptr = boost::make_shared<VectorDouble>();
 
  979    post_proc_volume->getOpPtrVector().push_back(
 
  981    post_proc_volume->getOpPtrVector().push_back(
 
  983    post_proc_volume->getOpPtrVector().push_back(
 
  988    post_proc_volume->getOpPtrVector().push_back(
 
  992            post_proc_volume->getPostProcMesh(),
 
  993            post_proc_volume->getMapGaussPts(),
 
  995            {{
"U", u_ptr}, {
"ERROR", e_ptr}},
 
  997            {{
"GRAD", grad_ptr}},
 
 
 1014      boost::function<
double(
const double, 
const double, 
const double)> f_u,
 
 1015      boost::function<
double(
const double, 
const double, 
const double)>
 
 1017      boost::function<
double(
const double)> 
a,
 
 1018      boost::function<
double(
const double)> diff_a,
 
 1019      boost::shared_ptr<ForcesAndSourcesCore> &domain_lhs_fe,
 
 1020      boost::shared_ptr<ForcesAndSourcesCore> &boundary_lhs_fe,
 
 1021      boost::shared_ptr<ForcesAndSourcesCore> &domain_rhs_fe,
 
 1022      boost::shared_ptr<ForcesAndSourcesCore> &boundary_rhs_fe,
 
 1025      bool trans = 
true)
 const {
 
 1029    domain_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
 
 1031    boundary_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
 
 1033    domain_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
 
 1035    boundary_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
 
 1039    domain_lhs_fe->getRuleHook = vol_rule;
 
 1040    domain_rhs_fe->getRuleHook = vol_rule;
 
 1041    boundary_lhs_fe->getRuleHook = face_rule;
 
 1042    boundary_rhs_fe->getRuleHook = face_rule;
 
 1047    boost::shared_ptr<VectorDouble> values_at_integration_ptr =
 
 1048        boost::make_shared<VectorDouble>();
 
 1050    boost::shared_ptr<MatrixDouble> grad_at_integration_ptr =
 
 1051        boost::make_shared<MatrixDouble>();
 
 1053    boost::shared_ptr<VectorDouble> multiplier_at_integration_ptr =
 
 1054        boost::make_shared<VectorDouble>();
 
 1058    domain_lhs_fe->getOpPtrVector().push_back(
 
 1061    domain_lhs_fe->getOpPtrVector().push_back(
 
 1064    domain_lhs_fe->getOpPtrVector().push_back(
new OpKt(
 
 1065        a, diff_a, values_at_integration_ptr, grad_at_integration_ptr));
 
 1068    domain_rhs_fe->getOpPtrVector().push_back(
 
 1071    domain_rhs_fe->getOpPtrVector().push_back(
 
 1074    domain_rhs_fe->getOpPtrVector().push_back(
new OpResF_Domain(
 
 1075        f_source, 
a, values_at_integration_ptr, grad_at_integration_ptr));
 
 1078    boundary_lhs_fe->getOpPtrVector().push_back(
new OpC(trans));
 
 1081    boundary_rhs_fe->getOpPtrVector().push_back(
 
 1084    boundary_rhs_fe->getOpPtrVector().push_back(
 
 1087    boundary_rhs_fe->getOpPtrVector().push_back(
 
 1088        new OpRes_g(f_u, values_at_integration_ptr));
 
 1089    boundary_rhs_fe->getOpPtrVector().push_back(
 
 
 
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#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 ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
const double c
speed of light (cm/ns)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
PostProcBrokenMeshInMoab< VolumeElementForcesAndSourcesCore > PostProcFE
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr auto field_name
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
const VectorDouble & getFieldData() const
Get DOF values on entity.
const VectorDofs & getFieldDofs() const
Get DOF data structures (const version)
const VectorInt & getIndices() const
Get global indices of degrees of freedom on entity.
default operator for TRI element
boost::function< int(int order_row, int order_col, int order_data)> RuleHookFun
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Specialization for double precision scalar field values calculation.
Operator for inverting matrices at integration points.
Post post-proc data at points from hash maps.
Set inverse jacobian to base functions.
Volume finite element base.
Create finite elements instances.
MoFEM::Interface & mField
MoFEMErrorCode creatFEToPostProcessResults(boost::shared_ptr< PostProcFE > &post_proc_volume) const
Create finite element to post-process results.
CreateFiniteElements(MoFEM::Interface &m_field)
MoFEMErrorCode createFEToAssembleMatrixAndVectorForNonlinearProblem(boost::function< double(const double, const double, const double)> f_u, boost::function< double(const double, const double, const double)> f_source, boost::function< double(const double)> a, boost::function< double(const double)> diff_a, boost::shared_ptr< ForcesAndSourcesCore > &domain_lhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &boundary_lhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &domain_rhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &boundary_rhs_fe, ForcesAndSourcesCore::RuleHookFun vol_rule, ForcesAndSourcesCore::RuleHookFun face_rule=FaceRule(), bool trans=true) const
Create finite element to calculate matrix and vectors.
MoFEMErrorCode createFEToEvaluateError(boost::function< double(const double, const double, const double)> f_u, boost::function< FTensor::Tensor1< double, 3 >(const double, const double, const double)> g_u, Vec global_error, boost::shared_ptr< ForcesAndSourcesCore > &domain_error) const
Create finite element to calculate error.
MoFEMErrorCode createFEToAssembleMatrixAndVector(boost::function< double(const double, const double, const double)> f_u, boost::function< double(const double, const double, const double)> f_source, boost::shared_ptr< ForcesAndSourcesCore > &domain_lhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &boundary_lhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &domain_rhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &boundary_rhs_fe, bool trans=true) const
Create finite element to calculate matrix and vectors.
Set integration rule to boundary elements.
int operator()(int p_row, int p_col, int p_data) const
template class for integration oh the right hand side
int nbIntegrationPts
number of integration points
virtual MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data)=0
Class dedicated to assemble operator to global system vector.
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
This function is called by finite element.
OpBaseRhs(const std::string field_name)
virtual MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)=0
Class dedicated to integrate operator.
Calculate constrains matrix.
MatrixDouble locMat
local constrains matrix
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
integrate local constrains matrix
OpC(const bool assemble_transpose)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate local constrains matrix.
int nbIntegrationPts
number of integration points
const bool assembleTranspose
assemble transpose, i.e. CT if set to true
int nbCols
number of dofs on column
Vec globalError
ghost vector with global (integrated over volume) error
UVal uValue
function with exact solution
OpError(UVal u_value, GVal g_value, boost::shared_ptr< VectorDouble > &field_vals, boost::shared_ptr< MatrixDouble > &grad_vals, Vec global_error)
boost::shared_ptr< MatrixDouble > gradVals
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
boost::function< double(const double, const double, const double)> UVal
FTensor::Index< 'i', 3 > i
GVal gValue
function with exact solution for gradient
boost::shared_ptr< VectorDouble > fieldVals
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate error.
boost::function< FTensor::Tensor1< double, 3 >(const double, const double, const double)> GVal
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data)
Assemble error.
Operator calculate source term,.
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data)
assemble local entity vector to the global right hand side
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate local entity vector.
boost::function< double(const double, const double, const double)> FSource
Calculate the grad-grad operator and assemble matrix.
virtual MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate grad-grad operator.
int nbCols
number if dof on column
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Do calculations for give operator.
bool isDiag
true if this block is on diagonal
MatrixDouble locMat
local entity block matrix
int nbIntegrationPts
number of integration points
FTensor::Index< 'i', 3 > i
summit Index
virtual MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Assemble local entity block matrix.
boost::function< double(const double)> diffA
boost::shared_ptr< MatrixDouble > gradVals
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate grad-grad operator.
boost::shared_ptr< VectorDouble > fieldVals
boost::function< double(const double)> A
OpKt(boost::function< double(const double)> a, boost::function< double(const double)> diff_a, boost::shared_ptr< VectorDouble > &field_vals, boost::shared_ptr< MatrixDouble > &grad_vals)
boost::shared_ptr< VectorDouble > lambdaVals
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate local constrains vector.
OpResF_Boundary(boost::shared_ptr< VectorDouble > &lambda_vals)
boost::shared_ptr< MatrixDouble > gradVals
boost::shared_ptr< VectorDouble > fieldVals
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate local entity vector.
boost::function< double(const double)> A
FTensor::Index< 'i', 3 > i
OpResF_Domain(FSource f_source, boost::function< double(const double)> a, boost::shared_ptr< VectorDouble > &field_vals, boost::shared_ptr< MatrixDouble > &grad_vals)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate local constrains vector.
boost::shared_ptr< VectorDouble > fieldVals
OpRes_g(FVal f_value, boost::shared_ptr< VectorDouble > &field_vals)
Assemble constrains vector.
FTensor::Number< 2 > NZ
z-direction index
Op_g(FVal f_value, const string field_name="L", const double beta=1)
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data)
assemble constrains vectors
boost::function< double(const double, const double, const double)> FVal
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate local constrains vector.
FTensor::Number< 0 > NX
x-direction index
FTensor::Number< 1 > NY
y-direction index
FVal fValue
Function pointer evaluating values of "U" at the boundary.
Set integration rule to volume elements.
int operator()(int, int, int p) const