7#ifndef __POISSONOPERATORS_HPP__
8#define __POISSONOPERATORS_HPP__
28struct OpK :
public VolumeElementForcesAndSourcesCore::UserDataOperator {
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;
305struct OpC :
public FaceElementForcesAndSourcesCore::UserDataOperator {
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();
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();
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();
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.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpK
const double c
speed of light (cm/ns)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
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 dofs values
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
const VectorInt & getIndices() const
Get global indices of dofs on entity.
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.
Get value at integration points for scalar field.
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)
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
boost::function< double(const double, const double, const double)> UVal
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
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.
boost::function< double(const double, const double, const double)> FVal
Set integration rule to volume elements.
int operator()(int, int, int p) const