v0.13.0
poisson_2d_nonhomogeneous.hpp
Go to the documentation of this file.
1 #ifndef __POISSON2DNONHOMOGENEOUS_HPP__
2 #define __POISSON2DNONHOMOGENEOUS_HPP__
3 
4 #include <stdlib.h>
5 #include <BasicFiniteElements.hpp>
6 
9 
12 
14 
16 
17 constexpr auto VecSetValues = MoFEM::VecSetValues<MoFEM::EssentialBcStorage>;
18 constexpr auto MatSetValues = MoFEM::MatSetValues<MoFEM::EssentialBcStorage>;
19 
21 
22 typedef boost::function<double(const double, const double, const double)>
24 
25 struct OpDomainLhs : public OpFaceEle {
26 public:
27  OpDomainLhs(std::string row_field_name, std::string col_field_name)
28  : OpFaceEle(row_field_name, col_field_name, OpFaceEle::OPROWCOL) {
29  sYmm = true;
30  }
31 
32  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
33  EntityType col_type, EntData &row_data,
34  EntData &col_data) {
36 
37  const int nb_row_dofs = row_data.getIndices().size();
38  const int nb_col_dofs = col_data.getIndices().size();
39 
40  if (nb_row_dofs && nb_col_dofs) {
41 
42  locLhs.resize(nb_row_dofs, nb_col_dofs, false);
43  locLhs.clear();
44 
45  // get element area
46  const double area = getMeasure();
47 
48  // get number of integration points
49  const int nb_integration_points = getGaussPts().size2();
50  // get integration weights
51  auto t_w = getFTensor0IntegrationWeight();
52 
53  // get derivatives of base functions on row
54  auto t_row_diff_base = row_data.getFTensor1DiffN<2>();
55 
56  // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL MATRIX
57  for (int gg = 0; gg != nb_integration_points; gg++) {
58 
59  const double a = t_w * area;
60 
61  for (int rr = 0; rr != nb_row_dofs; ++rr) {
62  // get derivatives of base functions on column
63  auto t_col_diff_base = col_data.getFTensor1DiffN<2>(gg, 0);
64 
65  for (int cc = 0; cc != nb_col_dofs; cc++) {
66 
67  locLhs(rr, cc) += t_row_diff_base(i) * t_col_diff_base(i) * a;
68 
69  // move to the derivatives of the next base functions on column
70  ++t_col_diff_base;
71  }
72 
73  // move to the derivatives of the next base functions on row
74  ++t_row_diff_base;
75  }
76 
77  // move to the weight of the next integration point
78  ++t_w;
79  }
80 
81  // FILL VALUES OF LOCAL MATRIX ENTRIES TO THE GLOBAL MATRIX
82 
83  // Fill value to local stiffness matrix ignoring boundary DOFs
84  CHKERR MatSetValues(getKSPB(), row_data, col_data, &locLhs(0, 0),
85  ADD_VALUES);
86 
87  // Fill values of symmetric local stiffness matrix
88  if (row_side != col_side || row_type != col_type) {
89  transLocLhs.resize(nb_col_dofs, nb_row_dofs, false);
90  noalias(transLocLhs) = trans(locLhs);
91  CHKERR MatSetValues(getKSPB(), col_data, row_data, &transLocLhs(0, 0),
92  ADD_VALUES);
93  }
94  }
95 
97  }
98 
99 private:
100  boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
102 };
103 
104 struct OpDomainRhs : public OpFaceEle {
105 public:
106  OpDomainRhs(std::string field_name, ScalarFunc source_term_function)
107  : OpFaceEle(field_name, OpFaceEle::OPROW),
108  sourceTermFunc(source_term_function) {}
109 
110  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
112 
113  const int nb_dofs = data.getIndices().size();
114 
115  if (nb_dofs) {
116 
117  locRhs.resize(nb_dofs, false);
118  locRhs.clear();
119 
120  // get element area
121  const double area = getMeasure();
122 
123  // get number of integration points
124  const int nb_integration_points = getGaussPts().size2();
125  // get integration weights
126  auto t_w = getFTensor0IntegrationWeight();
127  // get coordinates of the integration point
128  auto t_coords = getFTensor1CoordsAtGaussPts();
129 
130  // get base functions
131  auto t_base = data.getFTensor0N();
132 
133  // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL VECTOR
134  for (int gg = 0; gg != nb_integration_points; gg++) {
135 
136  const double a = t_w * area;
137  double body_source =
138  sourceTermFunc(t_coords(0), t_coords(1), t_coords(2));
139 
140  for (int rr = 0; rr != nb_dofs; rr++) {
141 
142  locRhs[rr] += t_base * body_source * a;
143 
144  // move to the next base function
145  ++t_base;
146  }
147 
148  // move to the weight of the next integration point
149  ++t_w;
150  // move to the coordinates of the next integration point
151  ++t_coords;
152  }
153 
154  // FILL VALUES OF THE GLOBAL VECTOR ENTRIES FROM THE LOCAL ONES
155  CHKERR VecSetValues(getKSPf(), data, &*locRhs.begin(), ADD_VALUES);
156  }
157 
159  }
160 
161 private:
164 };
165 
166 struct OpBoundaryLhs : public OpEdgeEle {
167 public:
168  OpBoundaryLhs(std::string row_field_name, std::string col_field_name)
169  : OpEdgeEle(row_field_name, col_field_name, OpEdgeEle::OPROWCOL) {
170  sYmm = true;
171  }
172 
173  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
174  EntityType col_type, EntData &row_data,
175  EntData &col_data) {
177 
178  const int nb_row_dofs = row_data.getIndices().size();
179  const int nb_col_dofs = col_data.getIndices().size();
180 
181  if (nb_row_dofs && nb_col_dofs) {
182 
183  locBoundaryLhs.resize(nb_row_dofs, nb_col_dofs, false);
184  locBoundaryLhs.clear();
185 
186  // get (boundary) element length
187  const double edge = getMeasure();
188 
189  // get number of integration points
190  const int nb_integration_points = getGaussPts().size2();
191  // get integration weights
192  auto t_w = getFTensor0IntegrationWeight();
193 
194  // get base functions on row
195  auto t_row_base = row_data.getFTensor0N();
196 
197  // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL MATRIX
198  for (int gg = 0; gg != nb_integration_points; gg++) {
199  const double a = t_w * edge;
200 
201  for (int rr = 0; rr != nb_row_dofs; ++rr) {
202  // get base functions on column
203  auto t_col_base = col_data.getFTensor0N(gg, 0);
204 
205  for (int cc = 0; cc != nb_col_dofs; cc++) {
206 
207  locBoundaryLhs(rr, cc) += t_row_base * t_col_base * a;
208 
209  // move to the next base functions on column
210  ++t_col_base;
211  }
212  // move to the next base functions on row
213  ++t_row_base;
214  }
215 
216  // move to the weight of the next integration point
217  ++t_w;
218  }
219 
220  // FILL VALUES OF LOCAL MATRIX ENTRIES TO THE GLOBAL MATRIX
221  CHKERR MatSetValues(getKSPB(), row_data, col_data, &locBoundaryLhs(0, 0),
222  ADD_VALUES);
223  if (row_side != col_side || row_type != col_type) {
224  transLocBoundaryLhs.resize(nb_col_dofs, nb_row_dofs, false);
225  noalias(transLocBoundaryLhs) = trans(locBoundaryLhs);
226  CHKERR MatSetValues(getKSPB(), col_data, row_data,
227  &transLocBoundaryLhs(0, 0), ADD_VALUES);
228  }
229  }
230 
232  }
233 
234 private:
236 };
237 
238 struct OpBoundaryRhs : public OpEdgeEle {
239 public:
240  OpBoundaryRhs(std::string field_name, ScalarFunc boundary_function)
241  : OpEdgeEle(field_name, OpEdgeEle::OPROW),
242  boundaryFunc(boundary_function) {}
243 
244  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
246 
247  const int nb_dofs = data.getIndices().size();
248 
249  if (nb_dofs) {
250 
251  locBoundaryRhs.resize(nb_dofs, false);
252  locBoundaryRhs.clear();
253 
254  // get (boundary) element length
255  const double edge = getMeasure();
256 
257  // get number of integration points
258  const int nb_integration_points = getGaussPts().size2();
259  // get integration weights
260  auto t_w = getFTensor0IntegrationWeight();
261  // get coordinates at integration point
262  auto t_coords = getFTensor1CoordsAtGaussPts();
263 
264  // get base function
265  auto t_base = data.getFTensor0N();
266 
267  // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL VECTOR
268  for (int gg = 0; gg != nb_integration_points; gg++) {
269 
270  const double a = t_w * edge;
271  double boundary_term =
272  boundaryFunc(t_coords(0), t_coords(1), t_coords(2));
273 
274  for (int rr = 0; rr != nb_dofs; rr++) {
275 
276  locBoundaryRhs[rr] += t_base * boundary_term * a;
277 
278  // move to the next base function
279  ++t_base;
280  }
281 
282  // move to the weight of the next integration point
283  ++t_w;
284  // move to the coordinates of the next integration point
285  ++t_coords;
286  }
287 
288  // FILL VALUES OF LOCAL VECTOR ENTRIES TO THE GLOBAL VECTOR
289  CHKERR VecSetValues(getKSPf(), data, &*locBoundaryRhs.begin(),
290  ADD_VALUES);
291  }
292 
294  }
295 
296 private:
299 };
300 
301 }; // namespace Poisson2DNonhomogeneousOperators
302 
303 #endif //__POISSON2DNONHOMOGENEOUS_HPP__
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
MoFEM::FaceElementForcesAndSourcesCore FaceEle
#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
boost::function< double(const double, const double, const double)> ScalarFunc
const double body_source
EntitiesFieldData::EntData EntData
bool sYmm
If true assume that matrix is symmetric structure.
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 VectorInt & getIndices() const
Get global indices of dofs on entity.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
auto getFTensor0IntegrationWeight()
Get integration weights.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
OpBoundaryLhs(std::string row_field_name, std::string col_field_name)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpBoundaryRhs(std::string field_name, ScalarFunc boundary_function)
OpDomainLhs(std::string row_field_name, std::string col_field_name)
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
OpDomainRhs(std::string field_name, ScalarFunc source_term_function)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.