v0.14.0
Loading...
Searching...
No Matches
poisson_2d_lagrange_multiplier.hpp
Go to the documentation of this file.
1#ifndef __POISSON2DLAGRANGEMULTIPLIER_HPP__
2#define __POISSON2DLAGRANGEMULTIPLIER_HPP__
3
4#include <stdlib.h>
6
9
12
14
16
18
19// const double body_source = 10;
20typedef boost::function<double(const double, const double, const double)>
22
23struct OpDomainLhsK : public OpFaceEle {
24public:
25 OpDomainLhsK(std::string row_field_name, std::string col_field_name,
26 boost::shared_ptr<std::vector<unsigned char>> boundary_marker = nullptr)
27 : OpFaceEle(row_field_name, col_field_name, OpFaceEle::OPROWCOL),
28 boundaryMarker(boundary_marker) {
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
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 CHKERR MatSetValues(getKSPB(), row_data, col_data, &locLhs(0, 0),
84 ADD_VALUES);
85
86 // Fill values of symmetric stiffness matrix in global system of equations
87 if (row_side != col_side || row_type != col_type) {
88 transLocLhs.resize(nb_col_dofs, nb_row_dofs, false);
89 noalias(transLocLhs) = trans(locLhs);
90 CHKERR MatSetValues(getKSPB(), col_data, row_data, &transLocLhs(0, 0),
91 ADD_VALUES);
92 }
93
94 }
95
97 }
98
99private:
100 boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
101 MatrixDouble locLhs, transLocLhs;
102};
103
104struct OpDomainRhsF : public OpFaceEle {
105public:
106 OpDomainRhsF(std::string field_name, ScalarFunc source_term_function,
107 boost::shared_ptr<std::vector<unsigned char>> boundary_marker = nullptr)
109 sourceTermFunc(source_term_function), boundaryMarker(boundary_marker) {}
110
111 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
113
114 const int nb_dofs = data.getIndices().size();
115
116 if (nb_dofs) {
117
118 locRhs.resize(nb_dofs, false);
119 locRhs.clear();
120
121 // get element area
122 const double area = getMeasure();
123
124 // get number of integration points
125 const int nb_integration_points = getGaussPts().size2();
126 // get integration weights
127 auto t_w = getFTensor0IntegrationWeight();
128 // get coordinates of the integration point
129 auto t_coords = getFTensor1CoordsAtGaussPts();
130
131 // get base functions
132 auto t_base = data.getFTensor0N();
133
134 // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL VECTOR
135 for (int gg = 0; gg != nb_integration_points; gg++) {
136
137 const double a = t_w * area;
138 double body_source =
139 sourceTermFunc(t_coords(0), t_coords(1), t_coords(2));
140
141 for (int rr = 0; rr != nb_dofs; rr++) {
142
143 locRhs[rr] += t_base * body_source * a;
144
145 // move to the next base function
146 ++t_base;
147 }
148
149 // move to the weight of the next integration point
150 ++t_w;
151 // move to the coordinates of the next integration point
152 ++t_coords;
153 }
154
155 // FILL VALUES OF THE GLOBAL VECTOR ENTRIES FROM THE LOCAL ONES
156
157 CHKERR VecSetValues(getKSPf(), data, &*locRhs.begin(), ADD_VALUES);
158 }
159
161 }
162
163private:
165 boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
166 VectorDouble locRhs;
167};
168
169struct OpBoundaryLhsC : public OpEdgeEle {
170public:
171 OpBoundaryLhsC(std::string row_field_name, std::string col_field_name)
172 : OpEdgeEle(row_field_name, col_field_name, OpEdgeEle::OPROWCOL) {
173 sYmm = false;
174 }
175
176 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
177 EntityType col_type, EntData &row_data,
178 EntData &col_data) {
180
181 const int nb_row_dofs = row_data.getIndices().size();
182 const int nb_col_dofs = col_data.getIndices().size();
183
184 if (nb_row_dofs && nb_col_dofs) {
185
186 locBoundaryLhs.resize(nb_row_dofs, nb_col_dofs, false);
187 locBoundaryLhs.clear();
188
189 // get (boundary) element length
190 const double edge = getMeasure();
191
192 // get number of integration points
193 const int nb_integration_points = getGaussPts().size2();
194 // get integration weights
195 auto t_w = getFTensor0IntegrationWeight();
196
197 // get base functions on row
198 auto t_row_base = row_data.getFTensor0N();
199
200 // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL MATRIX
201 for (int gg = 0; gg != nb_integration_points; gg++) {
202 const double a = t_w * edge;
203
204 for (int rr = 0; rr != nb_row_dofs; ++rr) {
205 // get base functions on column
206 auto t_col_base = col_data.getFTensor0N(gg, 0);
207
208 for (int cc = 0; cc != nb_col_dofs; cc++) {
209
210 locBoundaryLhs(rr, cc) += t_row_base * t_col_base * a;
211
212 // move to the next base functions on column
213 ++t_col_base;
214 }
215 // move to the next base functions on row
216 ++t_row_base;
217 }
218
219 // move to the weight of the next integration point
220 ++t_w;
221 }
222
223 // FILL VALUES OF LOCAL MATRIX ENTRIES TO THE GLOBAL MATRIX
224 CHKERR MatSetValues(getKSPB(), row_data, col_data, &locBoundaryLhs(0, 0),
225 ADD_VALUES);
226
227 }
228
230 }
231
232private:
234};
235
236struct OpBoundaryRhsG : public OpEdgeEle {
237public:
238 OpBoundaryRhsG(std::string field_name, ScalarFunc boundary_function)
240 boundaryFunc(boundary_function) {}
241
242 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
244
245 const int nb_dofs = data.getIndices().size();
246
247 if (nb_dofs) {
248
249 locBoundaryRhs.resize(nb_dofs, false);
250 locBoundaryRhs.clear();
251
252 // get (boundary) element length
253 const double edge = getMeasure();
254
255 // get number of integration points
256 const int nb_integration_points = getGaussPts().size2();
257 // get integration weights
258 auto t_w = getFTensor0IntegrationWeight();
259 // get coordinates at integration point
260 auto t_coords = getFTensor1CoordsAtGaussPts();
261
262 // get base function
263 auto t_base = data.getFTensor0N();
264
265 // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL VECTOR
266 for (int gg = 0; gg != nb_integration_points; gg++) {
267
268 const double a = t_w * edge;
269 double boundary_term =
270 boundaryFunc(t_coords(0), t_coords(1), t_coords(2));
271
272 for (int rr = 0; rr != nb_dofs; rr++) {
273
274 locBoundaryRhs[rr] += t_base * boundary_term * a;
275
276 // move to the next base function
277 ++t_base;
278 }
279
280 // move to the weight of the next integration point
281 ++t_w;
282 // move to the coordinates of the next integration point
283 ++t_coords;
284 }
285
286 // FILL VALUES OF LOCAL VECTOR ENTRIES TO THE GLOBAL VECTOR
287 CHKERR VecSetValues(getKSPf(), data, &*locBoundaryRhs.begin(),
288 ADD_VALUES);
289 }
290
292 }
293
294private:
296 VectorDouble locBoundaryRhs;
297};
298
299}; // namespace Poisson2DLagrangeMultiplierOperators
300
301#endif //__POISSON2DLAGRANGEMULTIPLIER_HPP__
constexpr double a
#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.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
boost::function< double(const double, const double, const double)> ScalarFunc
constexpr auto field_name
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.
double getMeasure() const
get measure of element
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
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)
Operator for bi-linear form, usually to calculate values on left hand side.
OpBoundaryLhsC(std::string row_field_name, std::string col_field_name)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpBoundaryRhsG(std::string field_name, ScalarFunc boundary_function)
OpDomainLhsK(std::string row_field_name, std::string col_field_name, boost::shared_ptr< std::vector< unsigned char > > boundary_marker=nullptr)
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.
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpDomainRhsF(std::string field_name, ScalarFunc source_term_function, boost::shared_ptr< std::vector< unsigned char > > boundary_marker=nullptr)