7 #ifndef __POISSONOPERATORS_HPP__
8 #define __POISSONOPERATORS_HPP__
12 using PostProcFE = PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore>;
47 nbRows = row_data.getIndices().size();
52 nbCols = col_data.getIndices().size();
59 if (row_side == col_side && row_type == col_type) {
96 double vol = getVolume();
98 auto t_w = getFTensor0IntegrationWeight();
100 auto t_row_grad = row_data.getFTensor1DiffN<3>();
104 const double alpha = t_w * vol;
110 for (
int rr = 0; rr !=
nbRows; rr++) {
112 auto t_col_grad = col_data.getFTensor1DiffN<3>(gg, 0);
114 for (
int cc = 0; cc !=
nbCols; cc++) {
116 a += alpha * (t_row_grad(
i) * t_col_grad(
i));
138 const int *row_indices = &*row_data.getIndices().data().begin();
140 const int *col_indices = &*col_data.getIndices().data().begin();
141 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
142 : getFEMethod()->snes_B;
145 &*
locMat.data().begin(), ADD_VALUES);
152 &*
locMat.data().begin(), ADD_VALUES);
176 nbRows = row_data.getIndices().size();
217 :
public OpBaseRhs<VolumeElementForcesAndSourcesCore::UserDataOperator> {
219 typedef boost::function<
double(
const double,
const double,
const double)>
246 double vol = getVolume();
248 auto t_w = getFTensor0IntegrationWeight();
250 auto t_v = data.getFTensor0N();
252 auto t_coords = getFTensor1CoordsAtGaussPts();
257 vol * t_w *
fSource(t_coords(
NX), t_coords(
NY), t_coords(
NZ));
262 for (
int rr = 0; rr !=
nbRows; rr++) {
282 const int *indices = &*data.getIndices().data().begin();
284 const double *vals = &*
locVec.data().begin();
285 Vec f = getFEMethod()->ksp_f != PETSC_NULL ? getFEMethod()->ksp_f
286 : getFEMethod()->snes_f;
304 OpC(
const bool assemble_transpose)
315 nbRows = row_data.getIndices().size();
320 nbCols = col_data.getIndices().size();
353 const double area = getArea();
355 auto t_w = getFTensor0IntegrationWeight();
357 auto t_row = row_data.getFTensor0N();
360 const double alpha = area * t_w;
365 for (
int rr = 0; rr !=
nbRows; rr++) {
367 auto t_col = col_data.getFTensor0N(gg, 0);
369 for (
int cc = 0; cc !=
nbCols; cc++) {
371 c += alpha * t_row * t_col;
389 const int *row_indices = &*row_data.getIndices().data().begin();
391 const int *col_indices = &*col_data.getIndices().data().begin();
392 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
393 : getFEMethod()->snes_B;
396 &*
locMat.data().begin(), ADD_VALUES);
402 &*
locMat.data().begin(), ADD_VALUES);
417 :
public OpBaseRhs<FaceElementForcesAndSourcesCore::UserDataOperator> {
419 typedef boost::function<
double(
const double,
const double,
const double)>
446 const double area = getArea() *
bEta;
448 auto t_w = getFTensor0IntegrationWeight();
450 auto t_l = data.getFTensor0N();
452 auto t_coords = getFTensor1CoordsAtGaussPts();
458 area * t_w *
fValue(t_coords(
NX), t_coords(
NY), t_coords(
NZ));
463 for (
int rr = 0; rr !=
nbRows; rr++) {
479 const int *indices = &*data.getIndices().data().begin();
480 const double *vals = &*
locVec.data().begin();
481 Vec f = getFEMethod()->ksp_f != PETSC_NULL ? getFEMethod()->ksp_f
482 : getFEMethod()->snes_f;
492 :
public OpBaseRhs<VolumeElementForcesAndSourcesCore::UserDataOperator> {
494 typedef boost::function<
double(
const double,
const double,
const double)>
496 typedef boost::function<FTensor::Tensor1<double, 3>(
497 const double,
const double,
const double)>
501 boost::shared_ptr<VectorDouble> &field_vals,
502 boost::shared_ptr<MatrixDouble> &grad_vals,
Vec global_error)
510 nbRows = row_data.getFieldData().size();
538 data.getFieldData().clear();
540 const double vol = getVolume();
542 auto t_w = getFTensor0IntegrationWeight();
546 auto t_grad = getFTensor1FromMat<3>(*
gradVals);
548 auto t_coords = getFTensor1CoordsAtGaussPts();
553 double alpha = vol * t_w;
555 double exact_u =
uValue(t_coords(
NX), t_coords(
NY), t_coords(
NZ));
557 t_exact_grad =
gValue(t_coords(
NX), t_coords(
NY), t_coords(
NZ));
559 t_error_grad(
i) = t_grad(
i) - t_exact_grad(
i);
561 double error = pow(t_u - exact_u, 2) + t_error_grad(
i) * t_error_grad(
i);
563 data.getFieldData()[0] += alpha * error;
578 data.getFieldDofs()[0]->getFieldData() = sqrt(data.getFieldData()[0]);
587 OpKt(boost::function<
double(
const double)>
a,
588 boost::function<
double(
const double)> diff_a,
589 boost::shared_ptr<VectorDouble> &field_vals,
590 boost::shared_ptr<MatrixDouble> &grad_vals)
609 double vol = getVolume();
611 auto t_w = getFTensor0IntegrationWeight();
615 auto t_grad = getFTensor1FromMat<3>(*
gradVals);
617 auto t_row_grad = row_data.getFTensor1DiffN<3>();
621 const double alpha = t_w * vol;
622 const double beta = alpha *
A(t_u);
624 t_gamma(
i) = (alpha *
diffA(t_u)) * t_grad(
i);
629 for (
int rr = 0; rr !=
nbRows; rr++) {
631 auto t_col = col_data.getFTensor0N(gg, 0);
633 auto t_col_grad = col_data.getFTensor1DiffN<3>(gg, 0);
635 for (
int cc = 0; cc !=
nbCols; cc++) {
637 a += (t_row_grad(
i) * beta) * t_col_grad(
i) +
638 t_row_grad(
i) * (t_gamma(
i) * t_col);
662 boost::shared_ptr<VectorDouble> &field_vals,
663 boost::shared_ptr<MatrixDouble> &grad_vals)
681 double vol = getVolume();
683 auto t_w = getFTensor0IntegrationWeight();
687 auto t_grad = getFTensor1FromMat<3>(*
gradVals);
689 auto t_v = data.getFTensor0N();
691 auto t_v_grad = data.getFTensor1DiffN<3>();
693 auto t_coords = getFTensor1CoordsAtGaussPts();
697 const double alpha = vol * t_w;
698 const double source_term =
701 grad_term(
i) = (alpha *
A(t_u)) * t_grad(
i);
706 for (
int rr = 0; rr !=
nbRows; rr++) {
708 t_a += t_v_grad(
i) * grad_term(
i) + t_v * source_term;
728 OpRes_g(
FVal f_value, boost::shared_ptr<VectorDouble> &field_vals)
744 const double area = getArea() *
bEta;
746 auto t_w = getFTensor0IntegrationWeight();
748 auto t_l = data.getFTensor0N();
752 auto t_coords = getFTensor1CoordsAtGaussPts();
757 double alpha = area * t_w;
761 for (
int rr = 0; rr !=
nbRows; rr++) {
763 (t_u -
fValue(t_coords(
NX), t_coords(
NY), t_coords(
NZ)));
793 const double area = getArea() *
bEta;
795 auto t_w = getFTensor0IntegrationWeight();
797 auto t_u = data.getFTensor0N();
804 double alpha = area * t_w;
808 for (
int rr = 0; rr !=
nbRows; rr++) {
809 t_a += alpha * t_u * t_lambda;
832 int operator()(
int,
int,
int p)
const {
return 2 * (p - 1); }
848 return 2 * p_data + 1;
866 boost::function<
double(
const double,
const double,
const double)> f_u,
867 boost::function<
double(
const double,
const double,
const double)>
869 boost::shared_ptr<ForcesAndSourcesCore> &domain_lhs_fe,
870 boost::shared_ptr<ForcesAndSourcesCore> &boundary_lhs_fe,
871 boost::shared_ptr<ForcesAndSourcesCore> &domain_rhs_fe,
872 boost::shared_ptr<ForcesAndSourcesCore> &boundary_rhs_fe,
873 bool trans =
true)
const {
877 domain_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
878 new VolumeElementForcesAndSourcesCore(
mField));
879 boundary_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
881 domain_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
882 new VolumeElementForcesAndSourcesCore(
mField));
883 boundary_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
887 domain_lhs_fe->getRuleHook =
VolRule();
888 domain_rhs_fe->getRuleHook =
VolRule();
889 boundary_lhs_fe->getRuleHook =
FaceRule();
890 boundary_rhs_fe->getRuleHook =
FaceRule();
894 domain_lhs_fe->getOpPtrVector().push_back(
new OpK());
896 domain_rhs_fe->getOpPtrVector().push_back(
new OpF(f_source));
898 boundary_lhs_fe->getOpPtrVector().push_back(
new OpC(trans));
900 boundary_rhs_fe->getOpPtrVector().push_back(
new Op_g(f_u));
909 boost::function<
double(
const double,
const double,
const double)> f_u,
914 boost::shared_ptr<ForcesAndSourcesCore> &domain_error)
const {
917 domain_error = boost::shared_ptr<ForcesAndSourcesCore>(
918 new VolumeElementForcesAndSourcesCore(
mField));
919 domain_error->getRuleHook =
VolRule();
923 boost::shared_ptr<VectorDouble> values_at_integration_ptr =
924 boost::make_shared<VectorDouble>();
926 boost::shared_ptr<MatrixDouble> grad_at_integration_ptr =
927 boost::make_shared<MatrixDouble>();
929 domain_error->getOpPtrVector().push_back(
930 new OpCalculateScalarFieldValues(
"U", values_at_integration_ptr));
932 domain_error->getOpPtrVector().push_back(
933 new OpCalculateScalarFieldGradient<3>(
"U", grad_at_integration_ptr));
935 domain_error->getOpPtrVector().push_back(
936 new OpError(f_u, g_u, values_at_integration_ptr,
937 grad_at_integration_ptr, global_error));
945 boost::shared_ptr<PostProcFE> &post_proc_volume)
const {
962 auto det_ptr = boost::make_shared<VectorDouble>();
963 auto jac_ptr = boost::make_shared<MatrixDouble>();
964 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
965 post_proc_volume->getOpPtrVector().push_back(
966 new OpCalculateHOJac<3>(jac_ptr));
967 post_proc_volume->getOpPtrVector().push_back(
968 new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
969 post_proc_volume->getOpPtrVector().push_back(
970 new OpSetHOInvJacToScalarBases<3>(
H1, inv_jac_ptr));
972 auto u_ptr = boost::make_shared<VectorDouble>();
973 auto grad_ptr = boost::make_shared<MatrixDouble>();
974 auto e_ptr = boost::make_shared<VectorDouble>();
976 post_proc_volume->getOpPtrVector().push_back(
977 new OpCalculateScalarFieldValues(
"U", u_ptr));
978 post_proc_volume->getOpPtrVector().push_back(
979 new OpCalculateScalarFieldGradient<3>(
"U", grad_ptr));
980 post_proc_volume->getOpPtrVector().push_back(
981 new OpCalculateScalarFieldValues(
"ERROR", e_ptr));
983 using OpPPMap = OpPostProcMapInMoab<3, 3>;
985 post_proc_volume->getOpPtrVector().push_back(
989 post_proc_volume->getPostProcMesh(),
990 post_proc_volume->getMapGaussPts(),
992 {{
"U", u_ptr}, {
"ERROR", e_ptr}},
994 {{
"GRAD", grad_ptr}},
1011 boost::function<
double(
const double,
const double,
const double)> f_u,
1012 boost::function<
double(
const double,
const double,
const double)>
1014 boost::function<
double(
const double)>
a,
1015 boost::function<
double(
const double)> diff_a,
1016 boost::shared_ptr<ForcesAndSourcesCore> &domain_lhs_fe,
1017 boost::shared_ptr<ForcesAndSourcesCore> &boundary_lhs_fe,
1018 boost::shared_ptr<ForcesAndSourcesCore> &domain_rhs_fe,
1019 boost::shared_ptr<ForcesAndSourcesCore> &boundary_rhs_fe,
1020 ForcesAndSourcesCore::RuleHookFun vol_rule,
1021 ForcesAndSourcesCore::RuleHookFun face_rule =
FaceRule(),
1022 bool trans =
true)
const {
1026 domain_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
1027 new VolumeElementForcesAndSourcesCore(mField));
1028 boundary_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
1030 domain_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
1031 new VolumeElementForcesAndSourcesCore(mField));
1032 boundary_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
1036 domain_lhs_fe->getRuleHook = vol_rule;
1037 domain_rhs_fe->getRuleHook = vol_rule;
1038 boundary_lhs_fe->getRuleHook = face_rule;
1039 boundary_rhs_fe->getRuleHook = face_rule;
1044 boost::shared_ptr<VectorDouble> values_at_integration_ptr =
1045 boost::make_shared<VectorDouble>();
1047 boost::shared_ptr<MatrixDouble> grad_at_integration_ptr =
1048 boost::make_shared<MatrixDouble>();
1050 boost::shared_ptr<VectorDouble> multiplier_at_integration_ptr =
1051 boost::make_shared<VectorDouble>();
1055 domain_lhs_fe->getOpPtrVector().push_back(
1056 new OpCalculateScalarFieldValues(
"U", values_at_integration_ptr));
1058 domain_lhs_fe->getOpPtrVector().push_back(
1059 new OpCalculateScalarFieldGradient<3>(
"U", grad_at_integration_ptr));
1061 domain_lhs_fe->getOpPtrVector().push_back(
new OpKt(
1062 a, diff_a, values_at_integration_ptr, grad_at_integration_ptr));
1065 domain_rhs_fe->getOpPtrVector().push_back(
1066 new OpCalculateScalarFieldValues(
"U", values_at_integration_ptr));
1068 domain_rhs_fe->getOpPtrVector().push_back(
1069 new OpCalculateScalarFieldGradient<3>(
"U", grad_at_integration_ptr));
1071 domain_rhs_fe->getOpPtrVector().push_back(
new OpResF_Domain(
1072 f_source,
a, values_at_integration_ptr, grad_at_integration_ptr));
1075 boundary_lhs_fe->getOpPtrVector().push_back(
new OpC(trans));
1078 boundary_rhs_fe->getOpPtrVector().push_back(
1079 new OpCalculateScalarFieldValues(
"U", values_at_integration_ptr));
1081 boundary_rhs_fe->getOpPtrVector().push_back(
1082 new OpCalculateScalarFieldValues(
"L", multiplier_at_integration_ptr));
1084 boundary_rhs_fe->getOpPtrVector().push_back(
1085 new OpRes_g(f_u, values_at_integration_ptr));
1086 boundary_rhs_fe->getOpPtrVector().push_back(
1098 #endif //__POISSONOPERATORS_HPP__