v0.13.0
PoissonOperators.hpp
Go to the documentation of this file.
1 #ifndef __POISSONOPERTATORS_HPP__
2 #define __POISSONOPERTATORS_HPP__
3 
4 #include <stdlib.h>
5 #include <BasicFiniteElements.hpp>
6 
7 namespace PoissonOps {
8 
11 
14 
16 
17 // Operators
18 
20 
21 const double body_source = 1;
22 
23 struct OpLhsUU : public OpVolEle {
24 public:
25  OpLhsUU(std::string domain_field1, std::string domain_field2)
26  : OpVolEle(domain_field1, domain_field2, OpVolEle::OPROWCOL) {
27  sYmm = true;
28  }
29 
30  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
31  EntityType col_type, EntData &row_data,
32  EntData &col_data) {
34 
35  const int nb_row_dofs = row_data.getIndices().size();
36  const int nb_col_dofs = col_data.getIndices().size();
37 
38  if (nb_row_dofs && nb_col_dofs) {
39  locLhs.resize(nb_row_dofs, nb_col_dofs, false);
40  locLhs.clear();
41  const int nb_integration_pts = getGaussPts().size2();
42  auto t_w = getFTensor0IntegrationWeight();
43 
44  auto t_row_diff_base = row_data.getFTensor1DiffN<2>();
45 
46  const double vol = getMeasure();
47  for (int gg = 0; gg != nb_integration_pts; ++gg) {
48  const double a = t_w * vol;
49  for (int rr = 0; rr != nb_row_dofs; ++rr) {
50  auto t_col_diff_base = col_data.getFTensor1DiffN<2>(gg, 0);
51  for (int cc = 0; cc != nb_col_dofs; ++cc) {
52  locLhs(rr, cc) += a * t_row_diff_base(i) * t_col_diff_base(i);
53  ++t_col_diff_base;
54  } // endFor cc
55  ++t_row_diff_base;
56  } // endFor rr
57  ++t_w;
58  } // endFor gg
59  CHKERR MatSetValues(getFEMethod()->ksp_B, row_data, col_data,
60  &locLhs(0, 0), ADD_VALUES);
61  if (row_side != col_side || row_type != col_type) {
62  transLocLhs.resize(nb_col_dofs, nb_row_dofs, false);
63  noalias(transLocLhs) = trans(locLhs);
64  CHKERR MatSetValues(getFEMethod()->ksp_B, col_data, row_data,
65  &transLocLhs(0, 0), ADD_VALUES);
66  } // (row_side != col_side || row_type != col_type)
67  } // endIf (nb_row_dofs && nb_col_dofs)
68 
70  }
71 
72 private:
74 };
75 
76 struct OpRhsU : public OpVolEle {
77 public:
78  OpRhsU(std::string domain_field) : OpVolEle(domain_field, OpVolEle::OPROW) {}
79  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
81 
82  const int nb_dofs = data.getIndices().size();
83  if (nb_dofs) {
84  locRhs.resize(nb_dofs, false);
85  locRhs.clear();
86  const int nb_integration_pts = getGaussPts().size2();
87 
88  auto t_base = data.getFTensor0N();
89  auto t_w = getFTensor0IntegrationWeight();
90  const double vol = getMeasure();
91  for (int gg = 0; gg != nb_integration_pts; ++gg) {
92  const double a = vol * t_w;
93 
94  for (int rr = 0; rr != nb_dofs; ++rr) {
95  locRhs[rr] += a * (t_base * body_source);
96  ++t_base;
97  } // endFor rr
98  ++t_w;
99  } // endFor gg
100  CHKERR VecSetOption(getFEMethod()->ksp_f, VEC_IGNORE_NEGATIVE_INDICES,
101  PETSC_TRUE);
102  CHKERR VecSetValues(getFEMethod()->ksp_f, data, &*locRhs.begin(),
103  ADD_VALUES);
104  } // endIf (nb_dofs)
105 
107  }
108 
109 private:
111 };
112 
113 }; // namespace PoissonOps
114 
115 #endif //__POISSONOPERTATORS_HPP__
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
EdgeElementForcesAndSourcesCoreSwitch< 0 > EdgeElementForcesAndSourcesCore
Edge finite element default.
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
FaceElementForcesAndSourcesCore::UserDataOperator OpVolEle
DataForcesAndSourcesCore::EntData EntData
FTensor::Index< 'i', 2 > i
const double body_source
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
OpLhsUU(std::string domain_field1, std::string domain_field2)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpRhsU(std::string domain_field)