v0.14.0
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 
9 namespace MoFEM {
10 
11 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, IntegrationType I,
12  typename OpBase>
14 
15 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, IntegrationType I,
16  typename OpBase>
18 
19 template <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 
34 protected:
35  boost::shared_ptr<MatrixDouble> yGradPtr;
37  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
38  EntitiesFieldData::EntData &col_data);
39 };
40 
41 template <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 
56 protected:
58  boost::shared_ptr<MatrixDouble> uPtr;
59  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
60  EntitiesFieldData::EntData &col_data);
61 };
62 
63 template <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 
78 protected:
80  boost::shared_ptr<MatrixDouble> yGradPtr;
81  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
82  EntitiesFieldData::EntData &col_data);
83 };
84 
85 template <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 
100 protected:
102  boost::shared_ptr<MatrixDouble> uPtr;
103  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
104  EntitiesFieldData::EntData &col_data);
105 };
106 
107 /**
108  * @brief Trilinear integrator form
109  * @ingroup mofem_forms
110  *
111  * @tparam EleOp
112  * @tparam A
113  * @tparam I
114  */
115 template <typename EleOp>
116 template <AssemblyType A>
117 template <IntegrationType I>
118 struct FormsIntegrators<EleOp>::Assembly<A>::TriLinearForm {
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>
129  using OpConvectiveTermLhsDu =
131 
132  template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM>
133  using OpConvectiveTermLhsDy =
135 };
136 
137 template <int SPACE_DIM, typename OpBase>
140  EntitiesFieldData::EntData &row_data,
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 
192 template <int SPACE_DIM, typename OpBase>
195  EntitiesFieldData::EntData &row_data,
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 
237 template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
240  EntitiesFieldData::EntData &row_data,
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 
296 template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
299  EntitiesFieldData::EntData &row_data,
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__
MoFEM::OpBaseImpl::nbIntegrationPts
int nbIntegrationPts
number of integration points
Definition: FormsIntegrators.hpp:238
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
MoFEM::OpConvectiveTermLhsDyImpl< 1, 1, SPACE_DIM, GAUSS, OpBase >::uPtr
boost::shared_ptr< MatrixDouble > uPtr
Definition: TriLinearFormsIntegrators.hpp:58
MoFEM::OpBaseImpl::locMat
MatrixDouble locMat
local entity block matrix
Definition: FormsIntegrators.hpp:249
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::ConstantFun
boost::function< double()> ConstantFun
Constant function type.
Definition: FormsIntegrators.hpp:166
OpBase
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
BASE_DIM
constexpr int BASE_DIM
Definition: dg_projection.cpp:15
MoFEM::OpConvectiveTermLhsDuImpl< 1, 1, SPACE_DIM, GAUSS, OpBase >::OpConvectiveTermLhsDuImpl
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;})
Definition: TriLinearFormsIntegrators.hpp:22
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
FIELD_DIM
constexpr int FIELD_DIM
Definition: child_and_parent.cpp:15
I
constexpr IntegrationType I
Definition: operators_tests.cpp:31
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1502
EleOp
OpBase
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:178
MoFEM::OpConvectiveTermLhsDyImpl< 1, 1, SPACE_DIM, GAUSS, OpBase >::OpConvectiveTermLhsDyImpl
OpConvectiveTermLhsDyImpl(const std::string field_name_row, const std::string field_name_col, boost::shared_ptr< MatrixDouble > u_ptr, ConstantFun alpha_fun=[]() { return 1;})
Definition: TriLinearFormsIntegrators.hpp:44
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
MoFEM::OpConvectiveTermLhsDuImpl< 1, FIELD_DIM, SPACE_DIM, GAUSS, OpBase >::alphaConstant
ConstantFun alphaConstant
Definition: TriLinearFormsIntegrators.hpp:79
MoFEM::OpConvectiveTermLhsDyImpl< 1, FIELD_DIM, SPACE_DIM, GAUSS, OpBase >::OpConvectiveTermLhsDyImpl
OpConvectiveTermLhsDyImpl(const std::string field_name_row, const std::string field_name_col, boost::shared_ptr< MatrixDouble > u_ptr, ConstantFun alpha_fun=[]() { return 1;})
Definition: TriLinearFormsIntegrators.hpp:88
MoFEM::OpConvectiveTermLhsDuImpl< 1, FIELD_DIM, SPACE_DIM, GAUSS, OpBase >::yGradPtr
boost::shared_ptr< MatrixDouble > yGradPtr
Definition: TriLinearFormsIntegrators.hpp:80
MoFEM::OpConvectiveTermLhsDuImpl< 1, 1, SPACE_DIM, GAUSS, OpBase >::yGradPtr
boost::shared_ptr< MatrixDouble > yGradPtr
Definition: TriLinearFormsIntegrators.hpp:35
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
MoFEM::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
MoFEM::OpConvectiveTermLhsDyImpl< 1, 1, SPACE_DIM, GAUSS, OpBase >::alphaConstant
ConstantFun alphaConstant
Definition: TriLinearFormsIntegrators.hpp:57
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::OpConvectiveTermLhsDuImpl
Definition: TriLinearFormsIntegrators.hpp:13
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
FTensor::Index< 'i', SPACE_DIM >
MoFEM::IntegrationType
IntegrationType
Form integrator integration types.
Definition: FormsIntegrators.hpp:136
MoFEM::OpBaseImpl::nbRowBaseFunctions
int nbRowBaseFunctions
number or row base functions
Definition: FormsIntegrators.hpp:239
MoFEM::OpBaseImpl::nbCols
int nbCols
number if dof on column
Definition: FormsIntegrators.hpp:237
MoFEM::OpConvectiveTermLhsDyImpl< 1, FIELD_DIM, SPACE_DIM, GAUSS, OpBase >::alphaConstant
ConstantFun alphaConstant
Definition: TriLinearFormsIntegrators.hpp:101
MoFEM::EntitiesFieldData::EntData::getFTensor1DiffN
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition: EntitiesFieldData.cpp:526
MoFEM::OpBaseImpl::nbRows
int nbRows
number of dofs on rows
Definition: FormsIntegrators.hpp:236
MoFEM::OpConvectiveTermLhsDuImpl< 1, FIELD_DIM, SPACE_DIM, GAUSS, OpBase >::OpConvectiveTermLhsDuImpl
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;})
Definition: TriLinearFormsIntegrators.hpp:66
MoFEM::OpMixTensorTimesGradImpl
Definition: BiLinearFormsIntegratorsImpl.hpp:386
MoFEM::OpConvectiveTermLhsDyImpl< 1, FIELD_DIM, SPACE_DIM, GAUSS, OpBase >::uPtr
boost::shared_ptr< MatrixDouble > uPtr
Definition: TriLinearFormsIntegrators.hpp:102
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::FormsIntegrators
Integrator forms.
Definition: FormsIntegrators.hpp:306
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEM::OpConvectiveTermLhsDuImpl< 1, 1, SPACE_DIM, GAUSS, OpBase >::alphaConstant
ConstantFun alphaConstant
Definition: TriLinearFormsIntegrators.hpp:36
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MoFEM::OpConvectiveTermLhsDyImpl
Definition: TriLinearFormsIntegrators.hpp:17