21 #ifndef __POISSONOPERATORS_HPP__
22 #define __POISSONOPERATORS_HPP__
59 nbRows = row_data.getIndices().size();
64 nbCols = col_data.getIndices().size();
71 if (row_side == col_side && row_type == col_type) {
109 double vol = getVolume();
111 auto t_w = getFTensor0IntegrationWeight();
113 auto t_row_grad = row_data.getFTensor1DiffN<3>();
117 const double alpha = t_w * vol;
123 for (
int rr = 0; rr !=
nbRows; rr++) {
125 auto t_col_grad = col_data.getFTensor1DiffN<3>(gg, 0);
127 for (
int cc = 0; cc !=
nbCols; cc++) {
129 a +=
alpha * (t_row_grad(
i) * t_col_grad(
i));
151 const int *row_indices = &*row_data.getIndices().data().begin();
153 const int *col_indices = &*col_data.getIndices().data().begin();
154 Mat
B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
155 : getFEMethod()->snes_B;
158 &*
locMat.data().begin(), ADD_VALUES);
165 &*
locMat.data().begin(), ADD_VALUES);
189 nbRows = row_data.getIndices().size();
230 :
public OpBaseRhs<VolumeElementForcesAndSourcesCore::UserDataOperator> {
232 typedef boost::function<double(
const double,
const double,
const double)>
259 double vol = getVolume();
261 auto t_w = getFTensor0IntegrationWeight();
263 auto t_v = data.getFTensor0N();
265 auto t_coords = getFTensor1CoordsAtGaussPts();
270 vol * t_w *
fSource(t_coords(
NX), t_coords(
NY), t_coords(
NZ));
275 for (
int rr = 0; rr !=
nbRows; rr++) {
295 const int *indices = &*data.getIndices().data().begin();
297 const double *vals = &*
locVec.data().begin();
298 Vec f = getFEMethod()->ksp_f != PETSC_NULL ? getFEMethod()->ksp_f
299 : getFEMethod()->snes_f;
317 OpC(
const bool assemble_transpose)
328 nbRows = row_data.getIndices().size();
333 nbCols = col_data.getIndices().size();
366 const double area = getArea();
368 auto t_w = getFTensor0IntegrationWeight();
370 auto t_row = row_data.getFTensor0N();
373 const double alpha = area * t_w;
378 for (
int rr = 0; rr !=
nbRows; rr++) {
380 auto t_col = col_data.getFTensor0N(gg, 0);
382 for (
int cc = 0; cc !=
nbCols; cc++) {
384 c +=
alpha * t_row * t_col;
402 const int *row_indices = &*row_data.getIndices().data().begin();
404 const int *col_indices = &*col_data.getIndices().data().begin();
405 Mat
B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
406 : getFEMethod()->snes_B;
409 &*
locMat.data().begin(), ADD_VALUES);
415 &*
locMat.data().begin(), ADD_VALUES);
430 :
public OpBaseRhs<FaceElementForcesAndSourcesCore::UserDataOperator> {
432 typedef boost::function<double(
const double,
const double,
const double)>
435 Op_g(
FVal f_value,
const string field_name =
"L",
const double beta = 1)
459 const double area = getArea() *
bEta;
461 auto t_w = getFTensor0IntegrationWeight();
463 auto t_l = data.getFTensor0N();
465 auto t_coords = getFTensor1CoordsAtGaussPts();
471 area * t_w *
fValue(t_coords(
NX), t_coords(
NY), t_coords(
NZ));
476 for (
int rr = 0; rr !=
nbRows; rr++) {
492 const int *indices = &*data.getIndices().data().begin();
493 const double *vals = &*
locVec.data().begin();
494 Vec f = getFEMethod()->ksp_f != PETSC_NULL ? getFEMethod()->ksp_f
495 : getFEMethod()->snes_f;
505 :
public OpBaseRhs<VolumeElementForcesAndSourcesCore::UserDataOperator> {
507 typedef boost::function<double(
const double,
const double,
const double)>
509 typedef boost::function<FTensor::Tensor1<double, 3>(
510 const double,
const double,
const double)>
514 boost::shared_ptr<VectorDouble> &field_vals,
515 boost::shared_ptr<MatrixDouble> &grad_vals,
Vec global_error)
523 nbRows = row_data.getFieldData().size();
551 data.getFieldData().clear();
553 const double vol = getVolume();
555 auto t_w = getFTensor0IntegrationWeight();
559 auto t_grad = getFTensor1FromMat<3>(*
gradVals);
561 auto t_coords = getFTensor1CoordsAtGaussPts();
566 double alpha = vol * t_w;
568 double exact_u =
uValue(t_coords(
NX), t_coords(
NY), t_coords(
NZ));
570 t_exact_grad =
gValue(t_coords(
NX), t_coords(
NY), t_coords(
NZ));
572 t_error_grad(
i) = t_grad(
i) - t_exact_grad(
i);
574 double error = pow(t_u - exact_u, 2) + t_error_grad(
i) * t_error_grad(
i);
576 data.getFieldData()[0] +=
alpha * error;
591 data.getFieldDofs()[0]->getFieldData() = sqrt(data.getFieldData()[0]);
600 OpKt(boost::function<
double(
const double)>
a,
601 boost::function<
double(
const double)> diff_a,
602 boost::shared_ptr<VectorDouble> &field_vals,
603 boost::shared_ptr<MatrixDouble> &grad_vals)
622 double vol = getVolume();
624 auto t_w = getFTensor0IntegrationWeight();
628 auto t_grad = getFTensor1FromMat<3>(*
gradVals);
630 auto t_row_grad = row_data.getFTensor1DiffN<3>();
634 const double alpha = t_w * vol;
635 const double beta =
alpha *
A(t_u);
642 for (
int rr = 0; rr !=
nbRows; rr++) {
644 auto t_col = col_data.getFTensor0N(gg, 0);
646 auto t_col_grad = col_data.getFTensor1DiffN<3>(gg, 0);
648 for (
int cc = 0; cc !=
nbCols; cc++) {
650 a += (t_row_grad(
i) * beta) * t_col_grad(
i) +
651 t_row_grad(
i) * (t_gamma(
i) * t_col);
666 boost::function<double(
const double)>
A;
667 boost::function<double(
const double)>
diffA;
675 boost::shared_ptr<VectorDouble> &field_vals,
676 boost::shared_ptr<MatrixDouble> &grad_vals)
694 double vol = getVolume();
696 auto t_w = getFTensor0IntegrationWeight();
700 auto t_grad = getFTensor1FromMat<3>(*
gradVals);
702 auto t_v = data.getFTensor0N();
704 auto t_v_grad = data.getFTensor1DiffN<3>();
706 auto t_coords = getFTensor1CoordsAtGaussPts();
710 const double alpha = vol * t_w;
711 const double source_term =
714 grad_term(
i) = (
alpha *
A(t_u)) * t_grad(
i);
719 for (
int rr = 0; rr !=
nbRows; rr++) {
721 t_a += t_v_grad(
i) * grad_term(
i) + t_v * source_term;
734 boost::function<double(
const double)>
A;
741 OpRes_g(
FVal f_value, boost::shared_ptr<VectorDouble> &field_vals)
757 const double area = getArea() *
bEta;
759 auto t_w = getFTensor0IntegrationWeight();
761 auto t_l = data.getFTensor0N();
765 auto t_coords = getFTensor1CoordsAtGaussPts();
770 double alpha = area * t_w;
774 for (
int rr = 0; rr !=
nbRows; rr++) {
776 (t_u -
fValue(t_coords(
NX), t_coords(
NY), t_coords(
NZ)));
806 const double area = getArea() *
bEta;
808 auto t_w = getFTensor0IntegrationWeight();
810 auto t_u = data.getFTensor0N();
817 double alpha = area * t_w;
821 for (
int rr = 0; rr !=
nbRows; rr++) {
822 t_a +=
alpha * t_u * t_lambda;
861 return 2 * p_data + 1;
879 boost::function<
double(
const double,
const double,
const double)> f_u,
880 boost::function<
double(
const double,
const double,
const double)>
882 boost::shared_ptr<ForcesAndSourcesCore> &domain_lhs_fe,
883 boost::shared_ptr<ForcesAndSourcesCore> &boundary_lhs_fe,
884 boost::shared_ptr<ForcesAndSourcesCore> &domain_rhs_fe,
885 boost::shared_ptr<ForcesAndSourcesCore> &boundary_rhs_fe,
886 bool trans =
true)
const {
890 domain_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
892 boundary_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
894 domain_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
896 boundary_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
900 domain_lhs_fe->getRuleHook =
VolRule();
901 domain_rhs_fe->getRuleHook =
VolRule();
902 boundary_lhs_fe->getRuleHook =
FaceRule();
903 boundary_rhs_fe->getRuleHook =
FaceRule();
907 domain_lhs_fe->getOpPtrVector().push_back(
new OpK());
909 domain_rhs_fe->getOpPtrVector().push_back(
new OpF(f_source));
911 boundary_lhs_fe->getOpPtrVector().push_back(
new OpC(trans));
913 boundary_rhs_fe->getOpPtrVector().push_back(
new Op_g(f_u));
922 boost::function<
double(
const double,
const double,
const double)> f_u,
927 boost::shared_ptr<ForcesAndSourcesCore> &domain_error)
const {
930 domain_error = boost::shared_ptr<ForcesAndSourcesCore>(
932 domain_error->getRuleHook =
VolRule();
936 boost::shared_ptr<VectorDouble> values_at_integration_ptr =
937 boost::make_shared<VectorDouble>();
939 boost::shared_ptr<MatrixDouble> grad_at_integration_ptr =
940 boost::make_shared<MatrixDouble>();
942 domain_error->getOpPtrVector().push_back(
943 new OpCalculateScalarFieldValues(
"U", values_at_integration_ptr));
945 domain_error->getOpPtrVector().push_back(
946 new OpCalculateScalarFieldGradient<3>(
"U", grad_at_integration_ptr));
948 domain_error->getOpPtrVector().push_back(
949 new OpError(f_u, g_u, values_at_integration_ptr,
950 grad_at_integration_ptr, global_error));
958 boost::shared_ptr<ForcesAndSourcesCore> &post_proc_volume)
const {
971 post_proc_volume = boost::shared_ptr<ForcesAndSourcesCore>(
974 CHKERR boost::static_pointer_cast<PostProcVolumeOnRefinedMesh>(
976 ->generateReferenceElementMesh();
977 CHKERR boost::static_pointer_cast<PostProcVolumeOnRefinedMesh>(
979 ->addFieldValuesPostProc(
"U");
980 CHKERR boost::static_pointer_cast<PostProcVolumeOnRefinedMesh>(
982 ->addFieldValuesPostProc(
"ERROR");
983 CHKERR boost::static_pointer_cast<PostProcVolumeOnRefinedMesh>(
985 ->addFieldValuesGradientPostProc(
"U");
994 boost::function<
double(
const double,
const double,
const double)> f_u,
995 boost::function<
double(
const double,
const double,
const double)>
997 boost::function<
double(
const double)>
a,
998 boost::function<
double(
const double)> diff_a,
999 boost::shared_ptr<ForcesAndSourcesCore> &domain_lhs_fe,
1000 boost::shared_ptr<ForcesAndSourcesCore> &boundary_lhs_fe,
1001 boost::shared_ptr<ForcesAndSourcesCore> &domain_rhs_fe,
1002 boost::shared_ptr<ForcesAndSourcesCore> &boundary_rhs_fe,
1003 ForcesAndSourcesCore::RuleHookFun vol_rule,
1004 ForcesAndSourcesCore::RuleHookFun face_rule =
FaceRule(),
1005 bool trans =
true)
const {
1009 domain_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
1011 boundary_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
1013 domain_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
1015 boundary_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
1019 domain_lhs_fe->getRuleHook = vol_rule;
1020 domain_rhs_fe->getRuleHook = vol_rule;
1021 boundary_lhs_fe->getRuleHook = face_rule;
1022 boundary_rhs_fe->getRuleHook = face_rule;
1027 boost::shared_ptr<VectorDouble> values_at_integration_ptr =
1028 boost::make_shared<VectorDouble>();
1030 boost::shared_ptr<MatrixDouble> grad_at_integration_ptr =
1031 boost::make_shared<MatrixDouble>();
1033 boost::shared_ptr<VectorDouble> multiplier_at_integration_ptr =
1034 boost::make_shared<VectorDouble>();
1038 domain_lhs_fe->getOpPtrVector().push_back(
1039 new OpCalculateScalarFieldValues(
"U", values_at_integration_ptr));
1041 domain_lhs_fe->getOpPtrVector().push_back(
1042 new OpCalculateScalarFieldGradient<3>(
"U", grad_at_integration_ptr));
1044 domain_lhs_fe->getOpPtrVector().push_back(
new OpKt(
1045 a, diff_a, values_at_integration_ptr, grad_at_integration_ptr));
1048 domain_rhs_fe->getOpPtrVector().push_back(
1049 new OpCalculateScalarFieldValues(
"U", values_at_integration_ptr));
1051 domain_rhs_fe->getOpPtrVector().push_back(
1052 new OpCalculateScalarFieldGradient<3>(
"U", grad_at_integration_ptr));
1054 domain_rhs_fe->getOpPtrVector().push_back(
new OpResF_Domain(
1055 f_source,
a, values_at_integration_ptr, grad_at_integration_ptr));
1058 boundary_lhs_fe->getOpPtrVector().push_back(
new OpC(trans));
1061 boundary_rhs_fe->getOpPtrVector().push_back(
1062 new OpCalculateScalarFieldValues(
"U", values_at_integration_ptr));
1064 boundary_rhs_fe->getOpPtrVector().push_back(
1065 new OpCalculateScalarFieldValues(
"L", multiplier_at_integration_ptr));
1067 boundary_rhs_fe->getOpPtrVector().push_back(
1068 new OpRes_g(f_u, values_at_integration_ptr));
1069 boundary_rhs_fe->getOpPtrVector().push_back(
ForcesAndSourcesCore::UserDataOperator UserDataOperator
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpK
#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.
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
UBlasVector< double > VectorDouble
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
constexpr auto VecSetValues
constexpr auto MatSetValues
DataForcesAndSourcesCore::EntData EntData
Deprecated interface functions.
Create finite elements instances.
MoFEM::Interface & mField
MoFEMErrorCode creatFEToPostProcessResults(boost::shared_ptr< ForcesAndSourcesCore > &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
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.
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)
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)
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)
Set integration rule to volume elements.
int operator()(int, int, int p) const