v0.15.0
Loading...
Searching...
No Matches
TriLinearFormsIntegrators.hpp
Go to the documentation of this file.
1/** \file TriLinearFormsIntegrators.hpp
2 * \brief Trilinear forms integrators
3 * \ingroup mofem_form
4 */
5
6#ifndef __TRILINEAR_FORMS_INTEGRATORS_HPP__
7#define __TRILINEAR_FORMS_INTEGRATORS_HPP__
8
9namespace MoFEM {
10
11template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, IntegrationType I,
12 typename OpBase>
14
15template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, IntegrationType I,
16 typename OpBase>
18
19template <int SPACE_DIM, typename OpBase>
21 : public OpBase {
23 const std::string field_name_row, const std::string field_name_col,
24 boost::shared_ptr<MatrixDouble> y_grad_ptr,
25 ConstantFun alpha_fun = []() { return 1; })
26 : OpBase(field_name_row, field_name_col, OpBase::OPROWCOL),
27 yGradPtr(y_grad_ptr), alphaConstant(alpha_fun) {
28
29 this->assembleTranspose = false;
30 this->onlyTranspose = false;
31 this->sYmm = false;
32 }
33
34protected:
35 boost::shared_ptr<MatrixDouble> yGradPtr;
39};
40
41template <int SPACE_DIM, typename OpBase>
43 : public OpBase {
45 const std::string field_name_row, const std::string field_name_col,
46 boost::shared_ptr<MatrixDouble> u_ptr,
47 ConstantFun alpha_fun = []() { return 1; })
48 : OpBase(field_name_row, field_name_col, OpBase::OPROWCOL), uPtr(u_ptr),
49 alphaConstant(alpha_fun) {
50
51 this->assembleTranspose = false;
52 this->onlyTranspose = false;
53 this->sYmm = false;
54 }
55
56protected:
58 boost::shared_ptr<MatrixDouble> uPtr;
61};
62
63template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
65 : public OpBase {
67 const std::string field_name_row, const std::string field_name_col,
68 boost::shared_ptr<MatrixDouble> y_grad_ptr,
69 ConstantFun alpha_fun = []() { return 1; })
70 : OpBase(field_name_row, field_name_col, OpBase::OPROWCOL),
71 yGradPtr(y_grad_ptr), alphaConstant(alpha_fun) {
72
73 this->assembleTranspose = false;
74 this->onlyTranspose = false;
75 this->sYmm = false;
76 }
77
78protected:
80 boost::shared_ptr<MatrixDouble> yGradPtr;
83};
84
85template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
87 : public OpBase {
89 const std::string field_name_row, const std::string field_name_col,
90 boost::shared_ptr<MatrixDouble> u_ptr,
91 ConstantFun alpha_fun = []() { return 1; })
92 : OpBase(field_name_row, field_name_col, OpBase::OPROWCOL), uPtr(u_ptr),
93 alphaConstant(alpha_fun) {
94
95 this->assembleTranspose = false;
96 this->onlyTranspose = false;
97 this->sYmm = false;
98 }
99
100protected:
102 boost::shared_ptr<MatrixDouble> uPtr;
103 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
105};
106
107/**
108 * @brief Trilinear integrator form
109 * @ingroup mofem_forms
110 *
111 * @tparam EleOp
112 * @tparam A
113 * @tparam I
114 */
115template <typename EleOp>
116template <AssemblyType A>
117template <IntegrationType I>
119
120 /**
121 * @brief Integrate \f$(\lambda_{ij},u_{i,j})_\Omega\f$
122 *
123 * @tparam SPACE_DIM
124 */
125 template <int SPACE_DIM>
127
128 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM>
131
132 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM>
135};
136
137template <int SPACE_DIM, typename OpBase>
141 EntitiesFieldData::EntData &col_data) {
143
144 // get element volume
145 const double vol = OpBase::getMeasure();
146 // get integration weights
147 auto t_w = OpBase::getFTensor0IntegrationWeight();
148 // get base function gradient on rows
149 auto t_row_base = row_data.getFTensor0N();
150
151 auto get_t_vec = [&](const int rr) {
152 std::array<double *, SPACE_DIM> ptrs;
153 for (auto i = 0; i != SPACE_DIM; ++i)
154 ptrs[i] = &OpBase::locMat(rr, i);
156 ptrs);
157 };
158
159 auto t_grad_y = getFTensor1FromMat<SPACE_DIM>(*yGradPtr);
161
162 const double alpha_constant = alphaConstant();
163 // loop over integration points
164 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
165 // take into account Jacobian
166 const double alpha = t_w * vol * alpha_constant;
167 // access local matrix
168 auto t_vec = get_t_vec(0);
169 // loop over rows base functions
170 int rr = 0;
171 for (; rr != OpBase::nbRows; rr++) {
172 // get column base functions gradient at gauss point gg
173 auto t_col_base = col_data.getFTensor0N(gg, 0);
174 // loop over columns
175 for (int cc = 0; cc != OpBase::nbCols / SPACE_DIM; cc++) {
176 // calculate element of local matrix
177 t_vec(i) += alpha * t_row_base * t_col_base * t_grad_y(i);
178 ++t_col_base;
179 ++t_vec;
180 }
181 ++t_row_base;
182 }
183 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
184 ++t_row_base;
185
186 ++t_grad_y;
187 ++t_w; // move to another integration weight
188 }
190};
191
192template <int SPACE_DIM, typename OpBase>
196 EntitiesFieldData::EntData &col_data) {
198
199 // get element volume
200 const double vol = OpBase::getMeasure();
201 // get integration weights
202 auto t_w = OpBase::getFTensor0IntegrationWeight();
203 // get base function gradient on rows
204 auto t_row_base = row_data.getFTensor0N();
205
206 auto t_u = getFTensor1FromMat<SPACE_DIM>(*uPtr);
208 const double alpha_constant = alphaConstant();
209 // loop over integration points
210 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
211 // take into account Jacobian
212 const double alpha = t_w * vol * alpha_constant;
213 // loop over rows base functions
214 auto a_mat_ptr = &*OpBase::locMat.data().begin();
215 int rr = 0;
216 for (; rr != OpBase::nbRows; rr++) {
217 // get column base functions gradient at gauss point gg
218 auto t_diff_col_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
219 // loop over columns
220 for (int cc = 0; cc != OpBase::nbCols; cc++) {
221 // calculate element of local matrix
222 (*a_mat_ptr) += alpha * t_row_base * t_diff_col_base(i) * t_u(i);
223 ++t_diff_col_base;
224 ++a_mat_ptr;
225 }
226 ++t_row_base;
227 }
228 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
229 ++t_row_base;
230
231 ++t_u;
232 ++t_w; // move to another integration weight
233 }
235};
236
237template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
241 EntitiesFieldData::EntData &col_data) {
243
244 // get element volume
245 const double vol = OpBase::getMeasure();
246 // get integration weights
247 auto t_w = OpBase::getFTensor0IntegrationWeight();
248 // get base function gradient on rows
249 auto t_row_base = row_data.getFTensor0N();
250
251 auto t_grad_y = getFTensor2FromMat<FIELD_DIM, SPACE_DIM>(*yGradPtr);
254
255 auto get_t_mat = [&](const int rr) {
256 std::array<double *, FIELD_DIM * SPACE_DIM> ptrs;
257 int s = 0;
258 for (int j = 0; j != FIELD_DIM; ++j)
259 for (auto i = 0; i != SPACE_DIM; ++i, ++s)
260 ptrs[s] = &OpBase::locMat(rr + j, i);
262 SPACE_DIM>(ptrs);
263 };
264
265 const double alpha_constant = alphaConstant();
266 // loop over integration points
267 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
268 // take into account Jacobian
269 const double alpha = t_w * vol * alpha_constant;
270
271 // loop over rows base functions
272 int rr = 0;
273 for (; rr != OpBase::nbRows / FIELD_DIM; rr++) {
274 // get column base functions gradient at gauss point gg
275 auto t_col_base = col_data.getFTensor0N(gg, 0);
276 // get mat
277 auto t_mat = get_t_mat(FIELD_DIM * rr);
278 // loop over columns
279 for (int cc = 0; cc != OpBase::nbCols / SPACE_DIM; cc++) {
280 // calculate element of local matrix
281 t_mat(I, k) += (alpha * t_row_base * t_col_base) * t_grad_y(I, k);
282 ++t_col_base;
283 ++t_mat;
284 }
285 ++t_row_base;
286 }
287 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
288 ++t_row_base;
289
290 ++t_grad_y;
291 ++t_w; // move to another integration weight
292 }
294};
295
296template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
300 EntitiesFieldData::EntData &col_data) {
302
303 // get element volume
304 const double vol = OpBase::getMeasure();
305 // get integration weights
306 auto t_w = OpBase::getFTensor0IntegrationWeight();
307 // get base function gradient on rows
308 auto t_row_base = row_data.getFTensor0N();
309
310 auto get_t_mat = [&](const int rr) {
311 std::array<double *, FIELD_DIM * FIELD_DIM> ptrs;
312 int s = 0;
313 for (int i = 0; i != FIELD_DIM; ++i)
314 for (int j = 0; j != FIELD_DIM; ++j, ++s)
315 ptrs[s] = &(OpBase::locMat(rr + i, j));
317 FIELD_DIM>(ptrs);
318 };
319
320 auto t_u = getFTensor1FromMat<SPACE_DIM>(*uPtr);
322
326
327 const double alpha_constant = alphaConstant();
328 // loop over integration points
329 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
330 // take into account Jacobian
331 const double alpha = t_w * vol * alpha_constant;
332
333 // loop over rows base functions
334 int rr = 0;
335 for (; rr != OpBase::nbRows / FIELD_DIM; ++rr) {
336 // get matrix vec
337 auto t_mat = get_t_mat(FIELD_DIM * rr);
338 // get column base functions gradient at gauss point gg
339 auto t_diff_col_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
340 // loop over columns
341 for (int cc = 0; cc != OpBase::nbCols / FIELD_DIM; ++cc) {
342 t_mat(I, L) +=
343 alpha * t_row_base * t_kd(I, L) * (t_diff_col_base(k) * t_u(k));
344 ++t_mat;
345 ++t_diff_col_base;
346 }
347 ++t_row_base;
348 }
349 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
350 ++t_row_base;
351
352 ++t_u;
353 ++t_w; // move to another integration weight
354 }
356};
357
358} // namespace MoFEM
359
360#endif //__TRILINEAR_FORMS_INTEGRATORS_HPP__
constexpr int SPACE_DIM
constexpr int FIELD_DIM
Kronecker Delta class symmetric.
#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()
constexpr int BASE_DIM
constexpr auto t_kd
IntegrationType
Form integrator integration types.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
boost::function< double()> ConstantFun
Constant function type.
constexpr IntegrationType I
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.
int nbRows
number of dofs on rows
int nbIntegrationPts
number of integration points
MatrixDouble locMat
local entity block matrix
int nbCols
number if dof on column
int nbRowBaseFunctions
number or row base functions
OpConvectiveTermLhsDuImpl(const std::string field_name_row, const std::string field_name_col, boost::shared_ptr< MatrixDouble > y_grad_ptr, ConstantFun alpha_fun=[]() { return 1;})
OpConvectiveTermLhsDuImpl(const std::string field_name_row, const std::string field_name_col, boost::shared_ptr< MatrixDouble > y_grad_ptr, ConstantFun alpha_fun=[]() { return 1;})
OpConvectiveTermLhsDyImpl(const std::string field_name_row, const std::string field_name_col, boost::shared_ptr< MatrixDouble > u_ptr, ConstantFun alpha_fun=[]() { return 1;})
OpConvectiveTermLhsDyImpl(const std::string field_name_row, const std::string field_name_col, boost::shared_ptr< MatrixDouble > u_ptr, ConstantFun alpha_fun=[]() { return 1;})