v0.14.0
BiLinearFormsIntegratorsImpl.hpp
Go to the documentation of this file.
1 /** \file BiLinearFormsIntegratorsImpl.hpp
2  * \brief Bilinear forms integrators (implementation)
3  * \ingroup mofem_form
4 
5  \todo Some operators could be optimised. To do that, we need to write tests
6  and use Valgrind to profile code, shaking cache misses. For example, some
7  operators should have iteration first over columns, then rows. ome operators.
8  Since those operators are used in many problems, an implementation must be
9  efficient.
10 
11 */
12 
13 #ifndef __BILINEAR_FORMS_INTEGRATORS_IMPL_HPP__
14 #define __BILINEAR_FORMS_INTEGRATORS_IMPL_HPP__
15 
16 namespace MoFEM {
17 
18 using CacheMatsTypeType =
19  std::map<std::pair<int, int>, boost::shared_ptr<MatrixDouble>>;
20 
21 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, IntegrationType I,
22  typename OpBase>
23 struct OpGradGradImpl {};
24 
25 template <int SPACE_DIM, typename OpBase>
26 struct OpGradGradImpl<1, 1, SPACE_DIM, GAUSS, OpBase> : public OpBase {
27  FTensor::Index<'i', SPACE_DIM> i; ///< summit Index
28  OpGradGradImpl(const std::string row_field_name,
29  const std::string col_field_name,
30  ScalarFun beta = scalar_fun_one,
31  boost::shared_ptr<Range> ents_ptr = nullptr)
32  : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
33  betaCoeff(beta) {
34  if (row_field_name == col_field_name)
35  this->sYmm = true;
36  }
37 
38 protected:
39  ScalarFun betaCoeff;
40  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
41  EntitiesFieldData::EntData &col_data);
42 };
43 
44 template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
45 struct OpGradGradImpl<1, FIELD_DIM, SPACE_DIM, GAUSS, OpBase> : public OpBase {
46  FTensor::Index<'i', SPACE_DIM> i; ///< summit Index
47  OpGradGradImpl(const std::string row_field_name,
48  const std::string col_field_name,
49  ScalarFun beta = scalar_fun_one,
50  boost::shared_ptr<Range> ents_ptr = nullptr)
51  : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
52  betaCoeff(beta) {
53  if (row_field_name == col_field_name)
54  this->sYmm = true;
55  }
56 
57 protected:
58  ScalarFun betaCoeff;
59  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
60  EntitiesFieldData::EntData &col_data);
61 };
62 
63 template <int BASE_DIM, int FIELD_DIM, IntegrationType I, typename OpBase>
64 struct OpMassImpl {};
65 
66 template <typename OpBase>
67 struct OpMassImpl<1, 1, GAUSS, OpBase> : public OpBase {
68 
69  OpMassImpl(
70  const std::string row_field_name, const std::string col_field_name,
71  ScalarFun beta = [](double, double, double) constexpr { return 1; },
72  boost::shared_ptr<Range> ents_ptr = nullptr,
73  boost::shared_ptr<MatrixDouble> cache_mat = nullptr)
74  : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
75  betaCoeff(beta) {
76  if (row_field_name == col_field_name)
77  this->sYmm = true;
78  }
79 
80 protected:
81  ScalarFun betaCoeff = [](double, double, double) constexpr { return 1; };
82  MoFEMErrorCode integrateImpl(EntitiesFieldData::EntData &row_data,
83  EntitiesFieldData::EntData &col_data,
84  double vol);
85  inline MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
86  EntitiesFieldData::EntData &col_data) {
87  return integrateImpl(row_data, col_data, OpBase::getMeasure());
88  }
89 };
90 
91 template <int FIELD_DIM, typename OpBase>
92 struct OpMassImpl<1, FIELD_DIM, GAUSS, OpBase>
93  : public OpMassImpl<1, 1, GAUSS, OpBase> {
94  using OpMassImpl<1, 1, GAUSS, OpBase>::OpMassImpl;
95 
96 protected:
97  using OpMassImpl<1, 1, GAUSS, OpBase>::betaCoeff;
98  MoFEMErrorCode integrateImpl(EntitiesFieldData::EntData &row_data,
99  EntitiesFieldData::EntData &col_data,
100  double vol);
101  inline MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
102  EntitiesFieldData::EntData &col_data) {
103  return integrateImpl(row_data, col_data, OpBase::getMeasure());
104  }
105 };
106 
107 template <int FIELD_DIM, typename OpBase>
108 struct OpMassImpl<3, FIELD_DIM, GAUSS, OpBase> : public OpBase {
109 
110  OpMassImpl(const std::string row_field_name, const std::string col_field_name,
111  ScalarFun beta = scalar_fun_one,
112  boost::shared_ptr<Range> ents_ptr = nullptr)
113  : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
114  betaCoeff(beta) {
115  if (row_field_name == col_field_name)
116  this->sYmm = true;
117  }
118 
119 protected:
120  ScalarFun betaCoeff;
121  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
122  EntitiesFieldData::EntData &col_data);
123 };
124 
125 template <typename OpBase>
126 struct OpMassImpl<3, 4, GAUSS, OpBase> : public OpBase {
127  OpMassImpl(const std::string row_field_name, const std::string col_field_name,
128  ScalarFun beta = scalar_fun_one,
129  boost::shared_ptr<Range> ents_ptr = nullptr)
130  : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
131  betaCoeff(beta) {
132  if (row_field_name == col_field_name)
133  this->sYmm = true;
134  }
135 
136 protected:
137  ScalarFun betaCoeff;
138  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
139  EntitiesFieldData::EntData &col_data);
140 };
141 
142 template <typename OpBase>
143 struct OpMassImpl<3, 9, GAUSS, OpBase> : public OpBase {
144  OpMassImpl(const std::string row_field_name, const std::string col_field_name,
145  ScalarFun beta = scalar_fun_one,
146  boost::shared_ptr<Range> ents_ptr = nullptr)
147  : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
148  betaCoeff(beta) {
149  if (row_field_name == col_field_name)
150  this->sYmm = true;
151  }
152 
153 protected:
154  ScalarFun betaCoeff;
155  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
156  EntitiesFieldData::EntData &col_data);
157 };
158 
159 template <int BASE_DIM, int FIELD_DIM, IntegrationType I, typename OpBase>
160 struct OpMassCacheImpl;
161 
162 template <int FIELD_DIM, IntegrationType I, typename OpBase>
163 struct OpMassCacheImpl<1, FIELD_DIM, I, OpBase>
164  : public OpMassImpl<1, FIELD_DIM, I, OpBase> {
165  OpMassCacheImpl(CacheMatsTypeType cache, const std::string row_field_name,
166  const std::string col_field_name, const double beta,
167  boost::shared_ptr<Range> ents_ptr = nullptr)
168  : OpMassImpl<1, FIELD_DIM, I, OpBase>(
169  row_field_name, col_field_name,
170  [](double, double, double) constexpr { return 1; }, ents_ptr),
171  cacheLocMats(cache), scalarBeta(beta) {}
172 
173  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
174  EntitiesFieldData::EntData &col_data);
175 
176 protected:
177  double scalarBeta;
178  CacheMatsTypeType cacheLocMats;
179 };
180 
181 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S, IntegrationType I,
182  typename OpBase>
183 struct OpGradSymTensorGradImpl {};
184 
185 template <int SPACE_DIM, int S, typename OpBase>
186 struct OpGradSymTensorGradImpl<1, SPACE_DIM, SPACE_DIM, S, GAUSS, OpBase>
187  : public OpBase {
188  FTensor::Index<'i', SPACE_DIM> i; ///< summit Index
189  OpGradSymTensorGradImpl(const std::string row_field_name,
190  const std::string col_field_name,
191  boost::shared_ptr<MatrixDouble> mat_D,
192  boost::shared_ptr<Range> ents_ptr = nullptr,
193  ScalarFun beta = scalar_fun_one)
194  : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
195  matD(mat_D), betaCoeff(beta) {
196  if (row_field_name == col_field_name)
197  this->sYmm = true;
198  }
199 
200 protected:
201  boost::shared_ptr<MatrixDouble> matD;
202  ScalarFun betaCoeff;
203  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
204  EntitiesFieldData::EntData &col_data);
205 };
206 
207 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S, IntegrationType I,
208  typename OpBase>
209 struct OpGradTensorGradImpl {};
210 
211 template <int SPACE_DIM, int S, typename OpBase>
212 struct OpGradTensorGradImpl<1, SPACE_DIM, SPACE_DIM, S, GAUSS, OpBase>
213  : public OpBase {
214  FTensor::Index<'i', SPACE_DIM> i; ///< summit Index
215  OpGradTensorGradImpl(const std::string row_field_name,
216  const std::string col_field_name,
217  boost::shared_ptr<MatrixDouble> mat_D,
218  ScalarFun beta = scalar_fun_one)
219  : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL), matD(mat_D),
220  betaCoeff(beta) {}
221 
222 protected:
223  boost::shared_ptr<MatrixDouble> matD;
224  ScalarFun betaCoeff;
225  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
226  EntitiesFieldData::EntData &col_data);
227 };
228 
229 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S, IntegrationType I,
230  typename OpBase>
231 struct OpGradGradSymTensorGradGradImpl {};
232 
233 template <int SPACE_DIM, int S, typename OpBase>
234 struct OpGradGradSymTensorGradGradImpl<1, 1, SPACE_DIM, S, GAUSS, OpBase>
235  : public OpBase {
236  FTensor::Index<'i', SPACE_DIM> i; ///< summit Index
237  OpGradGradSymTensorGradGradImpl(const std::string row_field_name,
238  const std::string col_field_name,
239  boost::shared_ptr<MatrixDouble> mat_D,
240  boost::shared_ptr<Range> ents_ptr = nullptr)
241  : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
242  matD(mat_D) {
243  if (row_field_name == col_field_name)
244  this->sYmm = true;
245  }
246 
247 protected:
248  boost::shared_ptr<MatrixDouble> matD;
249  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
250  EntitiesFieldData::EntData &col_data);
251 };
252 
253 template <int SPACE_DIM, IntegrationType I, typename OpBase>
254 struct OpMixDivTimesScalarImpl {};
255 
256 template <int SPACE_DIM, typename OpBase>
257 struct OpMixDivTimesScalarImpl<SPACE_DIM, GAUSS, OpBase> : public OpBase {
258  OpMixDivTimesScalarImpl(const std::string row_field_name,
259  const std::string col_field_name,
260  ConstantFun alpha_fun,
261  const bool assemble_transpose = false,
262  const bool only_transpose = false)
263  : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
264  alphaConstant(alpha_fun) {
265  this->assembleTranspose = assemble_transpose;
266  this->onlyTranspose = only_transpose;
267  }
268 
269 protected:
270  FTensor::Index<'i', SPACE_DIM> i; ///< summit Index
271  ConstantFun alphaConstant;
272  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
273  EntitiesFieldData::EntData &col_data);
274 };
275 
276 template <int SPACE_DIM, IntegrationType I, typename OpBase,
277  CoordinateTypes CoordSys>
278 struct OpMixDivTimesVecImpl {};
279 
280 template <int SPACE_DIM, typename OpBase, CoordinateTypes CoordSys>
281 struct OpMixDivTimesVecImpl<SPACE_DIM, GAUSS, OpBase, CoordSys>
282  : public OpBase {
283 
284  OpMixDivTimesVecImpl(const std::string row_field_name,
285  const std::string col_field_name, ConstantFun alpha_fun,
286  const bool assemble_transpose,
287  const bool only_transpose = false)
288  : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
289  alphaConstant(alpha_fun) {
290  this->assembleTranspose = assemble_transpose;
291  this->onlyTranspose = only_transpose;
292  }
293 
294  OpMixDivTimesVecImpl(const std::string row_field_name,
295  const std::string col_field_name, ConstantFun alpha_fun,
296  ScalarFun beta_fun, const bool assemble_transpose,
297  const bool only_transpose = false)
298  : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
299  alphaConstant(alpha_fun), betaCoeff(beta_fun) {
300  this->assembleTranspose = assemble_transpose;
301  this->onlyTranspose = only_transpose;
302  }
303 
304 protected:
305  FTensor::Index<'i', SPACE_DIM> i; ///< summit Index
306 
307  ConstantFun alphaConstant = []() constexpr { return 1; };
308  ScalarFun betaCoeff = [](double, double, double) constexpr { return 1; };
309 
310  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
311  EntitiesFieldData::EntData &col_data);
312 };
313 
314 template <int SPACE_DIM, IntegrationType I, typename OpBase,
315  CoordinateTypes COORDINATE_SYSTEM>
316 struct OpMixScalarTimesDivImpl {};
317 
318 template <int SPACE_DIM, typename OpBase, CoordinateTypes COORDINATE_SYSTEM>
319 struct OpMixScalarTimesDivImpl<SPACE_DIM, GAUSS, OpBase, COORDINATE_SYSTEM>
320  : public OpBase {
321  OpMixScalarTimesDivImpl(
322  const std::string row_field_name, const std::string col_field_name,
323  ScalarFun alpha_fun = [](double, double, double) constexpr { return 1; },
324  const bool assemble_transpose = false, const bool only_transpose = false)
325  : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
326  alphaConstant(alpha_fun) {
327  this->assembleTranspose = assemble_transpose;
328  this->onlyTranspose = only_transpose;
329  this->sYmm = false;
330  }
331 
332 protected:
333  FTensor::Index<'i', SPACE_DIM> i; ///< summit Index
334  ScalarFun alphaConstant;
335  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
336  EntitiesFieldData::EntData &col_data);
337 };
338 
339 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, IntegrationType I,
340  typename OpBase>
341 struct OpMixVectorTimesGradImpl;
342 
343 template <int SPACE_DIM, typename OpBase>
344 struct OpMixVectorTimesGradImpl<3, SPACE_DIM, SPACE_DIM, GAUSS, OpBase>
345  : public OpBase {
346  OpMixVectorTimesGradImpl(const std::string row_field_name,
347  const std::string col_field_name,
348  ConstantFun alpha_fun,
349  const bool assemble_transpose = false,
350  const bool only_transpose = false)
351  : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
352  alphaConstant(alpha_fun) {
353  this->assembleTranspose = assemble_transpose;
354  this->onlyTranspose = only_transpose;
355  }
356 
357 protected:
358  FTensor::Index<'i', SPACE_DIM> i; ///< summit Index
359  ConstantFun alphaConstant;
360  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
361  EntitiesFieldData::EntData &col_data);
362 };
363 
364 template <int SPACE_DIM, typename OpBase>
365 struct OpMixVectorTimesGradImpl<1, SPACE_DIM, SPACE_DIM, GAUSS, OpBase>
366  : public OpBase {
367  OpMixVectorTimesGradImpl(const std::string row_field_name,
368  const std::string col_field_name,
369  ConstantFun alpha_fun,
370  const bool assemble_transpose = false,
371  const bool only_transpose = false)
372  : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
373  alphaConstant(alpha_fun) {
374  this->assembleTranspose = assemble_transpose;
375  this->onlyTranspose = only_transpose;
376  }
377 
378 protected:
379  FTensor::Index<'i', SPACE_DIM> i; ///< summit Index
380  ConstantFun alphaConstant;
381  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
382  EntitiesFieldData::EntData &col_data);
383 };
384 
385 template <int SPACE_DIM, IntegrationType I, typename OpBase>
386 struct OpMixTensorTimesGradImpl {};
387 
388 template <int SPACE_DIM, typename OpBase>
389 struct OpMixTensorTimesGradImpl<SPACE_DIM, GAUSS, OpBase> : public OpBase {
390 
391  OpMixTensorTimesGradImpl(const std::string row_field_name,
392  const std::string col_field_name,
393  ConstantFun alpha_fun, const bool assemble_transpose,
394  const bool only_transpose = false)
395  : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
396  alphaConstant(alpha_fun) {
397  this->assembleTranspose = assemble_transpose;
398  this->onlyTranspose = only_transpose;
399  }
400 
401  OpMixTensorTimesGradImpl(const std::string row_field_name,
402  const std::string col_field_name,
403  ConstantFun alpha_fun, ScalarFun beta_coeff,
404  const bool assemble_transpose,
405  const bool only_transpose = false)
406  : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
407  alphaConstant(alpha_fun), betaCoeff(beta_coeff) {
408  this->assembleTranspose = assemble_transpose;
409  this->onlyTranspose = only_transpose;
410  }
411 
412 protected:
413  FTensor::Index<'i', SPACE_DIM> i; ///< summit Index
414  FTensor::Index<'j', SPACE_DIM> j; ///< summit Index
415  ConstantFun alphaConstant = []() constexpr { return 1; };
416  ScalarFun betaCoeff = [](double, double, double) constexpr { return 1; };
417  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
418  EntitiesFieldData::EntData &col_data);
419 };
420 
421 template <int SPACE_DIM, typename OpBase>
422 MoFEMErrorCode OpGradGradImpl<1, 1, SPACE_DIM, GAUSS, OpBase>::iNtegrate(
423  EntitiesFieldData::EntData &row_data,
424  EntitiesFieldData::EntData &col_data) {
425  MoFEMFunctionBegin;
426 
427  // get element volume
428  const double vol = OpBase::getMeasure();
429  // get integration weights
430  auto t_w = OpBase::getFTensor0IntegrationWeight();
431  // get base function gradient on rows
432  auto t_row_grad = row_data.getFTensor1DiffN<SPACE_DIM>();
433  // get coordinate at integration points
434  auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
435 
436  // loop over integration points
437  for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
438  const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
439  // take into account Jacobian
440  const double alpha = t_w * beta;
441  // loop over ros base functions
442  int rr = 0;
443  for (; rr != OpBase::nbRows; rr++) {
444  // get column base functions gradient at gauss point gg
445  auto t_col_grad = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
446  // loop over columns
447  for (int cc = 0; cc != OpBase::nbCols; cc++) {
448  // calculate element of local matrix
449  OpBase::locMat(rr, cc) += alpha * (t_row_grad(i) * t_col_grad(i));
450  ++t_col_grad; // move to another gradient of base function
451  // on column
452  }
453  ++t_row_grad; // move to another element of gradient of base
454  // function on row
455  }
456  for (; rr < OpBase::nbRowBaseFunctions; ++rr)
457  ++t_row_grad;
458 
459  ++t_coords;
460  ++t_w; // move to another integration weight
461  }
462  MoFEMFunctionReturn(0);
463 }
464 
465 template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
466 MoFEMErrorCode
467 OpGradGradImpl<1, FIELD_DIM, SPACE_DIM, GAUSS, OpBase>::iNtegrate(
468  EntitiesFieldData::EntData &row_data,
469  EntitiesFieldData::EntData &col_data) {
470  MoFEMFunctionBegin;
471 
472  // get element volume
473  const double vol = OpBase::getMeasure();
474  // get integration weights
475  auto t_w = OpBase::getFTensor0IntegrationWeight();
476  // get base function gradient on rows
477  auto t_row_grad = row_data.getFTensor1DiffN<SPACE_DIM>();
478  // get coordinate at integration points
479  auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
480 
481  FTensor::Index<'i', SPACE_DIM> i;
482  FTensor::Index<'j', FIELD_DIM> j;
483 
484  auto get_t_vec = [&](const int rr) {
485  std::array<double *, FIELD_DIM> ptrs;
486  for (auto i = 0; i != FIELD_DIM; ++i)
487  ptrs[i] = &OpBase::locMat(rr + i, i);
488  return FTensor::Tensor1<FTensor::PackPtr<double *, FIELD_DIM>, FIELD_DIM>(
489  ptrs);
490  };
491 
492  // loop over integration points
493  for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
494  const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
495  // take into account Jacobian
496  const double alpha = t_w * beta;
497  // loop over ros base functions
498  int rr = 0;
499  for (; rr != OpBase::nbRows / FIELD_DIM; rr++) {
500  // get diag vec
501  auto t_vec = get_t_vec(rr * FIELD_DIM);
502  // get column base functions gradient at gauss point gg
503  auto t_col_grad = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
504  // loop over columns
505  for (int cc = 0; cc != OpBase::nbCols / FIELD_DIM; cc++) {
506  // calculate element of local matrix
507  t_vec(j) += alpha * (t_row_grad(i) * t_col_grad(i));
508  ++t_col_grad; // move to another gradient of base function
509  // on column
510  ++t_vec;
511  }
512  ++t_row_grad; // move to another element of gradient of base
513  // function on row
514  }
515  for (; rr < OpBase::nbRowBaseFunctions; ++rr)
516  ++t_row_grad;
517 
518  ++t_coords;
519  ++t_w; // move to another integration weight
520  }
521  MoFEMFunctionReturn(0);
522 }
523 
524 template <typename OpBase>
525 MoFEMErrorCode OpMassImpl<1, 1, GAUSS, OpBase>::integrateImpl(
526  EntitiesFieldData::EntData &row_data,
527  EntitiesFieldData::EntData &col_data, double vol) {
528  MoFEMFunctionBegin;
529 
530 #ifndef NDEBUG
531  auto log_error = [&]() {
532  MOFEM_LOG("SELF", Sev::error) << "Row side " << OpBase::rowSide << " "
533  << CN::EntityTypeName(OpBase::rowType);
534  MOFEM_LOG("SELF", Sev::error) << "Col side " << OpBase::colSide << " "
535  << CN::EntityTypeName(OpBase::colType);
536  };
537 
538  if (row_data.getN().size2() < OpBase::nbRows) {
539  log_error();
540  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
541  "Wrong number of base functions on rows %d < %d",
542  row_data.getN().size2(), OpBase::nbRows);
543  }
544  if (col_data.getN().size2() < OpBase::nbCols) {
545  log_error();
546  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
547  "Wrong number of base functions on cols %d < %d",
548  col_data.getN().size2(), OpBase::nbCols);
549  }
550  if (row_data.getN().size1() != OpBase::nbIntegrationPts) {
551  log_error();
552  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
553  "Wrong number of integration points on rows %d != %d",
554  row_data.getN().size1(), OpBase::nbIntegrationPts);
555  }
556  if (col_data.getN().size1() != OpBase::nbIntegrationPts) {
557  log_error();
558  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
559  "Wrong number of integration points on cols %d != %d",
560  col_data.getN().size1(), OpBase::nbIntegrationPts);
561  }
562 #endif
563 
564  // get integration weights
565  auto t_w = OpBase::getFTensor0IntegrationWeight();
566  // get base function gradient on rows
567  auto t_row_base = row_data.getFTensor0N();
568  // get coordinate at integration points
569  auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
570  // loop over integration points
571  for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
572  const double beta = betaCoeff(t_coords(0), t_coords(1), t_coords(2));
573  // take into account Jacobian
574  const double alpha = t_w * beta;
575  // loop over rows base functions
576  auto a_mat_ptr = &*OpBase::locMat.data().begin();
577  int rr = 0;
578  for (; rr != OpBase::nbRows; rr++) {
579  // get column base functions gradient at gauss point gg
580  auto t_col_base = col_data.getFTensor0N(gg, 0);
581  // loop over columns
582  for (int cc = 0; cc != OpBase::nbCols; cc++) {
583  // calculate element of local matrix
584  *a_mat_ptr += alpha * (t_row_base * t_col_base);
585  ++t_col_base;
586  ++a_mat_ptr;
587  }
588  ++t_row_base;
589  }
590  for (; rr < OpBase::nbRowBaseFunctions; ++rr)
591  ++t_row_base;
592  ++t_coords;
593  ++t_w; // move to another integration weight
594  }
595 
596  OpBase::locMat *= vol;
597 
598  MoFEMFunctionReturn(0);
599 };
600 
601 template <int FIELD_DIM, typename OpBase>
602 MoFEMErrorCode OpMassImpl<1, FIELD_DIM, GAUSS, OpBase>::integrateImpl(
603  EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data,
604  double vol) {
605  MoFEMFunctionBegin;
606 
607  auto integrate = [&](auto &mat) {
608  MoFEMFunctionBeginHot;
609  // get integration weights
610  auto t_w = OpBase::getFTensor0IntegrationWeight();
611  // get base function gradient on rows
612  auto t_row_base = row_data.getFTensor0N();
613  // get coordinate at integration points
614  auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
615 
616  FTensor::Index<'i', FIELD_DIM> i;
617  auto get_t_vec = [&](const int rr) {
618  std::array<double *, FIELD_DIM> ptrs;
619  for (auto i = 0; i != FIELD_DIM; ++i)
620  ptrs[i] = &mat(rr + i, i);
621  return FTensor::Tensor1<FTensor::PackPtr<double *, FIELD_DIM>, FIELD_DIM>(
622  ptrs);
623  };
624 
625  // loop over integration points
626  for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
627  const double beta = betaCoeff(t_coords(0), t_coords(1), t_coords(2));
628  // take into account Jacobian
629  const double alpha = t_w * beta;
630  // loop over rows base functions
631  int rr = 0;
632  for (; rr != OpBase::nbRows / FIELD_DIM; rr++) {
633  // get column base functions gradient at gauss point gg
634  auto t_col_base = col_data.getFTensor0N(gg, 0);
635  // get mat vec
636  auto t_vec = get_t_vec(FIELD_DIM * rr);
637  // loop over columns
638  for (int cc = 0; cc != OpBase::nbCols / FIELD_DIM; cc++) {
639  // calculate element of local matrix
640  t_vec(i) += alpha * (t_row_base * t_col_base);
641  ++t_col_base;
642  ++t_vec;
643  }
644  ++t_row_base;
645  }
646  for (; rr < OpBase::nbRowBaseFunctions; ++rr)
647  ++t_row_base;
648  ++t_coords;
649  ++t_w; // move to another integration weight
650  }
651  MoFEMFunctionReturnHot(0);
652  };
653 
654  CHKERR integrate(OpBase::locMat);
655  OpBase::locMat *= vol;
656 
657  MoFEMFunctionReturn(0);
658 };
659 
660 template <int FIELD_DIM, typename OpBase>
661 MoFEMErrorCode OpMassImpl<3, FIELD_DIM, GAUSS, OpBase>::iNtegrate(
662  EntitiesFieldData::EntData &row_data,
663  EntitiesFieldData::EntData &col_data) {
664  FTensor::Index<'i', FIELD_DIM> i;
665  MoFEMFunctionBegin;
666  size_t nb_base_functions = row_data.getN().size2() / 3;
667  // // get element volume
668  const double vol = OpBase::getMeasure();
669  // get integration weights
670  auto t_w = OpBase::getFTensor0IntegrationWeight();
671  // get base function gradient on rows
672  auto t_row_base = row_data.getFTensor1N<3>();
673  // get coordinate at integration points
674  auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
675  // loop over integration points
676  for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
677  const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
678  // take into account Jacobian
679  const double alpha = t_w * beta;
680  // loop over rows base functions
681  auto a_mat_ptr = &*OpBase::locMat.data().begin();
682  int rr = 0;
683  for (; rr != OpBase::nbRows; rr++) {
684  // get column base functions gradient at gauss point gg
685  auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
686  // loop over columns
687  for (int cc = 0; cc != OpBase::nbCols; cc++) {
688  // calculate element of local matrix
689  (*a_mat_ptr) += alpha * (t_row_base(i) * t_col_base(i));
690  ++t_col_base;
691  ++a_mat_ptr;
692  }
693  ++t_row_base;
694  }
695  for (; rr < nb_base_functions; ++rr)
696  ++t_row_base;
697  ++t_coords;
698  ++t_w; // move to another integration weight
699  }
700  MoFEMFunctionReturn(0);
701 };
702 
703 template <typename OpBase>
704 MoFEMErrorCode OpMassImpl<3, 4, GAUSS, OpBase>::iNtegrate(
705  EntitiesFieldData::EntData &row_data,
706  EntitiesFieldData::EntData &col_data) {
707  MoFEMFunctionBegin;
708  FTensor::Index<'i', 2> I;
709  FTensor::Index<'k', 3> k;
710  auto get_t_vec = [&](const int rr) {
711  return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>{
712  &OpBase::locMat(rr + 0, 0), &OpBase::locMat(rr + 1, 1)};
713  };
714  size_t nb_base_functions = row_data.getN().size2() / 3;
715  // // get element volume
716  const double vol = OpBase::getMeasure();
717  // get integration weights
718  auto t_w = OpBase::getFTensor0IntegrationWeight();
719  // get base function gradient on rows
720  auto t_row_base = row_data.getFTensor1N<3>();
721  // get coordinate at integration points
722  auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
723  // loop over integration points
724  for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
725  const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
726  // take into account Jacobian
727  const double alpha = t_w * beta;
728  // loop over rows base functions
729  int rr = 0;
730  for (; rr != OpBase::nbRows / 2; rr++) {
731  // get column base functions gradient at gauss point gg
732  auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
733  auto t_vec = get_t_vec(2 * rr);
734  // loop over columns
735  for (int cc = 0; cc != OpBase::nbCols / 2; cc++) {
736  // calculate element of local matrix
737  t_vec(I) += alpha * (t_row_base(k) * t_col_base(k));
738  ++t_col_base;
739  ++t_vec;
740  }
741  ++t_row_base;
742  }
743  for (; rr < nb_base_functions; ++rr)
744  ++t_row_base;
745  ++t_coords;
746  ++t_w; // move to another integration weight
747  }
748  MoFEMFunctionReturn(0);
749 }
750 
751 template <typename OpBase>
752 MoFEMErrorCode OpMassImpl<3, 9, GAUSS, OpBase>::iNtegrate(
753  EntitiesFieldData::EntData &row_data,
754  EntitiesFieldData::EntData &col_data) {
755  MoFEMFunctionBegin;
756  FTensor::Index<'i', 3> i;
757  FTensor::Index<'k', 3> k;
758  auto get_t_vec = [&](const int rr) {
759  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>{
760  &OpBase::locMat(rr + 0, 0), &OpBase::locMat(rr + 1, 1),
761  &OpBase::locMat(rr + 2, 2)};
762  };
763  size_t nb_base_functions = row_data.getN().size2() / 3;
764  // // get element volume
765  const double vol = OpBase::getMeasure();
766  // get integration weights
767  auto t_w = OpBase::getFTensor0IntegrationWeight();
768  // get base function gradient on rows
769  auto t_row_base = row_data.getFTensor1N<3>();
770  // get coordinate at integration points
771  auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
772  // loop over integration points
773  for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
774  const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
775  // take into account Jacobian
776  const double alpha = t_w * beta;
777  // loop over rows base functions
778  int rr = 0;
779  for (; rr != OpBase::nbRows / 3; rr++) {
780  // get column base functions gradient at gauss point gg
781  auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
782  auto t_vec = get_t_vec(3 * rr);
783  // loop over columns
784  for (int cc = 0; cc != OpBase::nbCols / 3; cc++) {
785  // calculate element of local matrix
786  t_vec(i) += alpha * (t_row_base(k) * t_col_base(k));
787  ++t_col_base;
788  ++t_vec;
789  }
790  ++t_row_base;
791  }
792  for (; rr < nb_base_functions; ++rr)
793  ++t_row_base;
794  ++t_coords;
795  ++t_w; // move to another integration weight
796  }
797  MoFEMFunctionReturn(0);
798 };
799 
800 template <int FIELD_DIM, IntegrationType I, typename OpBase>
801 MoFEMErrorCode OpMassCacheImpl<1, FIELD_DIM, I, OpBase>::iNtegrate(
802  EntitiesFieldData::EntData &row_data,
803  EntitiesFieldData::EntData &col_data) {
804  MoFEMFunctionBegin;
805 
806  const auto vol = this->getMeasure();
807  const auto row_type = this->rowType;
808  const auto col_type = this->colType;
809  auto &loc_mat = this->locMat;
810 
811  auto p = std::make_pair(row_type, col_type);
812 
813  if (cacheLocMats[p]) {
814  if (cacheLocMats[p]->size1() != loc_mat.size1() &&
815  cacheLocMats[p]->size2() != loc_mat.size2()) {
816  cacheLocMats[p]->resize(loc_mat.size1(), loc_mat.size2());
817  CHKERR this->integrateImpl(row_data, col_data, 1);
818  *(cacheLocMats[p]) = loc_mat;
819  } else {
820  loc_mat = *(cacheLocMats[p]);
821  }
822  loc_mat *= scalarBeta * this->getMeasure();
823  } else {
824  CHKERR this->integrateImpl(row_data, col_data,
825  scalarBeta * this->getMeasure());
826  }
827  MoFEMFunctionReturn(0);
828 }
829 
830 template <int SPACE_DIM, int S, typename OpBase>
831 MoFEMErrorCode
832 OpGradSymTensorGradImpl<1, SPACE_DIM, SPACE_DIM, S, GAUSS, OpBase>::iNtegrate(
833  EntitiesFieldData::EntData &row_data,
834  EntitiesFieldData::EntData &col_data) {
835  MoFEMFunctionBegin;
836 
837  const size_t nb_row_dofs = row_data.getIndices().size();
838  const size_t nb_col_dofs = col_data.getIndices().size();
839 
840  if (nb_row_dofs && nb_col_dofs) {
841 
842  const bool diag = (row_data.getFieldEntities()[0]->getLocalUniqueId() ==
843  col_data.getFieldEntities()[0]->getLocalUniqueId());
844 
845  FTensor::Index<'i', SPACE_DIM> i;
846  FTensor::Index<'j', SPACE_DIM> j;
847  FTensor::Index<'k', SPACE_DIM> k;
848  FTensor::Index<'l', SPACE_DIM> l;
849 
850  // get element volume
851  double vol = OpBase::getMeasure();
852 
853  // get intergrayion weights
854  auto t_w = OpBase::getFTensor0IntegrationWeight();
855 
856  // get derivatives of base functions on rows
857  auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
858 
859  // Elastic stiffness tensor (4th rank tensor with minor and major
860  // symmetry)
861  auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, S>(*matD);
862 
863  // get coordinate at integration points
864  auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
865 
866  // iterate over integration points
867  for (int gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
868 
869  // calculate scalar weight times element volume
870  double a = t_w * vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
871 
872  // iterate over row base functions
873  int rr = 0;
874  for (; rr != OpBase::nbRows / SPACE_DIM; ++rr) {
875 
876  // get sub matrix for the row
877  auto t_m = OpBase::template getLocMat<SPACE_DIM>(SPACE_DIM * rr);
878 
879  FTensor::Christof<double, SPACE_DIM, SPACE_DIM> t_rowD;
880  // I mix up the indices here so that it behaves like a
881  // Dg. That way I don't have to have a separate wrapper
882  // class Christof_Expr, which simplifies things.
883  t_rowD(l, j, k) = t_D(i, j, k, l) * (a * t_row_diff_base(i));
884 
885  // get derivatives of base functions for columns
886  auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
887 
888  // iterate column base functions
889  const auto nb_cols = (diag) ? rr : OpBase::nbCols / SPACE_DIM - 1;
890  for (int cc = 0; cc <= nb_cols; ++cc) {
891 
892  // integrate block local stiffens matrix
893  t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base(k);
894 
895  // move to next column base function
896  ++t_col_diff_base;
897 
898  // move to next block of local stiffens matrix
899  ++t_m;
900  }
901 
902  // move to next row base function
903  ++t_row_diff_base;
904  }
905 
906  for (; rr < OpBase::nbRowBaseFunctions; ++rr)
907  ++t_row_diff_base;
908 
909  // move to next integration weight
910  ++t_w;
911  ++t_D;
912  ++t_coords;
913  }
914 
915  // Copy symmetry
916  if (diag) {
917  for (int rr = 0; rr != OpBase::nbRows / SPACE_DIM - 1; ++rr) {
918  auto t_m_rr = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
919  this->locMat, SPACE_DIM * rr, SPACE_DIM * (rr + 1));
920  auto t_m_cc = getFTensor2FromArray<SPACE_DIM, SPACE_DIM>(
921  this->locMat, SPACE_DIM * (rr + 1), SPACE_DIM * rr,
922  SPACE_DIM * OpBase::nbCols);
923  for (int cc = rr + 1; cc < OpBase::nbCols / SPACE_DIM; ++cc) {
924  t_m_rr(i, j) = t_m_cc(j, i);
925  ++t_m_rr;
926  ++t_m_cc;
927  }
928  }
929  }
930  }
931 
932  MoFEMFunctionReturn(0);
933 }
934 
935 template <int SPACE_DIM, int S, typename OpBase>
936 MoFEMErrorCode
937 OpGradGradSymTensorGradGradImpl<1, 1, SPACE_DIM, S, GAUSS, OpBase>::iNtegrate(
938  EntitiesFieldData::EntData &row_data,
939  EntitiesFieldData::EntData &col_data) {
940  MoFEMFunctionBegin;
941 
942  if (OpBase::nbRows && OpBase::nbCols) {
943 
944  FTensor::Index<'i', SPACE_DIM> i;
945  FTensor::Index<'j', SPACE_DIM> j;
946  FTensor::Index<'k', SPACE_DIM> k;
947  FTensor::Index<'l', SPACE_DIM> l;
948 
949  auto &row_hessian = row_data.getN(BaseDerivatives::SecondDerivative);
950  auto &col_hessian = col_data.getN(BaseDerivatives::SecondDerivative);
951 
952 #ifndef NDEBUG
953  if (row_hessian.size1() != OpBase::nbIntegrationPts) {
954  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
955  "Wrong number of integration pts (%d != %d)",
956  row_hessian.size1(), OpBase::nbIntegrationPts);
957  }
958  if (row_hessian.size2() !=
959  OpBase::nbRowBaseFunctions * SPACE_DIM * SPACE_DIM) {
960  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
961  "Wrong number of base functions (%d != %d)",
962  row_hessian.size2() / (SPACE_DIM * SPACE_DIM),
963  OpBase::nbRowBaseFunctions);
964  }
965  if (row_hessian.size2() < OpBase::nbRows * SPACE_DIM * SPACE_DIM) {
966  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
967  "Wrong number of base functions (%d < %d)", row_hessian.size2(),
968  OpBase::nbRows * SPACE_DIM * SPACE_DIM);
969  }
970  if (col_hessian.size1() != OpBase::nbIntegrationPts) {
971  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
972  "Wrong number of integration pts (%d != %d)",
973  col_hessian.size1(), OpBase::nbIntegrationPts);
974  }
975  if (col_hessian.size2() < OpBase::nbCols * SPACE_DIM * SPACE_DIM) {
976  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
977  "Wrong number of base functions (%d < %d)", col_hessian.size2(),
978  OpBase::nbRows * SPACE_DIM * SPACE_DIM);
979  }
980 #endif
981 
982  // get element volume
983  double vol = OpBase::getMeasure();
984 
985  // get intergrayion weights
986  auto t_w = OpBase::getFTensor0IntegrationWeight();
987 
988  auto t_row_diff2 = getFTensor2SymmetricLowerFromPtr<SPACE_DIM>(
989  &*row_hessian.data().begin());
990 
991  // Elastic stiffness tensor (4th rank tensor with minor and major
992  // symmetry)
993  auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, S>(*matD);
994 
995  // iterate over integration points
996  for (int gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
997 
998  // calculate scalar weight times element volume
999  double a = t_w * vol;
1000 
1001  // get sub matrix for the row
1002  auto m_ptr = &*OpBase::locMat.data().begin();
1003 
1004  // iterate over row base functions
1005  int rr = 0;
1006  for (; rr != OpBase::nbRows; ++rr) {
1007 
1008  FTensor::Tensor2_symmetric<double, SPACE_DIM> t_rowD;
1009  t_rowD(k, l) = t_D(i, j, k, l) * (a * t_row_diff2(i, j));
1010 
1011  // get derivatives of base functions for columns
1012  auto t_col_diff2 =
1013  getFTensor2SymmetricLowerFromPtr<SPACE_DIM>(&col_hessian(gg, 0));
1014 
1015  // iterate column base functions
1016  for (int cc = 0; cc != OpBase::nbCols; ++cc) {
1017 
1018  // integrate block local stiffens matrix
1019  *m_ptr += t_rowD(i, j) * t_col_diff2(i, j);
1020 
1021  // move to next column base function
1022  ++t_col_diff2;
1023 
1024  // move to next block of local stiffens matrix
1025  ++m_ptr;
1026  }
1027 
1028  // move to next row base function
1029  ++t_row_diff2;
1030  }
1031 
1032  for (; rr < OpBase::nbRowBaseFunctions; ++rr)
1033  ++t_row_diff2;
1034 
1035  // move to next integration weight
1036  ++t_w;
1037  ++t_D;
1038  }
1039  }
1040 
1041  MoFEMFunctionReturn(0);
1042 }
1043 
1044 template <int SPACE_DIM, int S, typename OpBase>
1045 MoFEMErrorCode
1046 OpGradTensorGradImpl<1, SPACE_DIM, SPACE_DIM, S, GAUSS, OpBase>::iNtegrate(
1047  EntitiesFieldData::EntData &row_data,
1048  EntitiesFieldData::EntData &col_data) {
1049  MoFEMFunctionBegin;
1050 
1051  const size_t nb_row_dofs = row_data.getIndices().size();
1052  const size_t nb_col_dofs = col_data.getIndices().size();
1053 
1054  if (nb_row_dofs && nb_col_dofs) {
1055 
1056  FTensor::Index<'i', SPACE_DIM> i;
1057  FTensor::Index<'j', SPACE_DIM> j;
1058  FTensor::Index<'k', SPACE_DIM> k;
1059  FTensor::Index<'l', SPACE_DIM> l;
1060 
1061  // get element volume
1062  double vol = OpBase::getMeasure();
1063 
1064  // get intergrayion weights
1065  auto t_w = OpBase::getFTensor0IntegrationWeight();
1066 
1067  // get derivatives of base functions on rows
1068  auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
1069 
1070  // stiffness tensor (4th rank tensor)
1071  auto t_D =
1072  getFTensor4FromMat<SPACE_DIM, SPACE_DIM, SPACE_DIM, SPACE_DIM, S>(
1073  *matD);
1074 
1075  // get coordinate at integration points
1076  auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
1077 
1078  // iterate over integration points
1079  for (int gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1080 
1081  // calculate scalar weight times element volume
1082  double a = t_w * vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
1083 
1084  // iterate over row base functions
1085  int rr = 0;
1086  for (; rr != OpBase::nbRows / SPACE_DIM; ++rr) {
1087 
1088  // get sub matrix for the row
1089  auto t_m = OpBase::template getLocMat<SPACE_DIM>(SPACE_DIM * rr);
1090 
1091  // calculate row
1092  FTensor::Tensor3<double, 3, 3, 3> t_row;
1093  t_row(i, k, l) = t_D(i, j, k, l) * (a * t_row_diff_base(j));
1094 
1095  // get derivatives of base functions for columns
1096  auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1097 
1098  // iterate column base functions
1099  for (int cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1100 
1101  // integrate block local stiffens matrix
1102  t_m(i, k) += t_row(i, k, l) * t_col_diff_base(l);
1103 
1104  // move to next column base function
1105  ++t_col_diff_base;
1106 
1107  // move to next block of local stiffens matrix
1108  ++t_m;
1109  }
1110 
1111  // move to next row base function
1112  ++t_row_diff_base;
1113  }
1114 
1115  for (; rr < OpBase::nbRowBaseFunctions; ++rr)
1116  ++t_row_diff_base;
1117 
1118  // move to next integration weight
1119  ++t_w;
1120  ++t_D;
1121  ++t_coords;
1122  }
1123  }
1124 
1125  MoFEMFunctionReturn(0);
1126 }
1127 
1128 template <int SPACE_DIM, typename OpBase>
1129 MoFEMErrorCode OpMixDivTimesScalarImpl<SPACE_DIM, GAUSS, OpBase>::iNtegrate(
1130  EntitiesFieldData::EntData &row_data,
1131  EntitiesFieldData::EntData &col_data) {
1132  MoFEMFunctionBegin;
1133 
1134  auto t_w = this->getFTensor0IntegrationWeight();
1135 
1136  size_t nb_base_functions = row_data.getN().size2() / 3;
1137  auto t_row_diff_base = row_data.getFTensor2DiffN<3, SPACE_DIM>();
1138 
1139  const double alpha_constant = alphaConstant();
1140 
1141  for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1142  const double alpha = alpha_constant * this->getMeasure() * t_w;
1143 
1144  size_t rr = 0;
1145  for (; rr != OpBase::nbRows; ++rr) {
1146  const double t_row_div_base = t_row_diff_base(i, i);
1147  auto t_col_base = col_data.getFTensor0N(gg, 0);
1148  for (size_t cc = 0; cc != OpBase::nbCols; ++cc) {
1149  this->locMat(rr, cc) += alpha * t_row_div_base * t_col_base;
1150  ++t_col_base;
1151  }
1152  ++t_row_diff_base;
1153  }
1154  for (; rr < nb_base_functions; ++rr)
1155  ++t_row_diff_base;
1156 
1157  ++t_w;
1158  }
1159 
1160  MoFEMFunctionReturn(0);
1161 }
1162 
1163 template <int SPACE_DIM, typename OpBase, CoordinateTypes CoordSys>
1164 MoFEMErrorCode
1165 OpMixDivTimesVecImpl<SPACE_DIM, GAUSS, OpBase, CoordSys>::iNtegrate(
1166  EntitiesFieldData::EntData &row_data,
1167  EntitiesFieldData::EntData &col_data) {
1168  MoFEMFunctionBegin;
1169 
1170  auto t_w = this->getFTensor0IntegrationWeight();
1171  auto t_coords = this->getFTensor1CoordsAtGaussPts();
1172 
1173  size_t nb_base_functions = row_data.getN().size2() / 3;
1174  auto t_row_diff_base = row_data.getFTensor2DiffN<3, SPACE_DIM>();
1175  auto t_row_base = row_data.getFTensor1N<3>();
1176  const double alpha_constant = alphaConstant() * this->getMeasure();
1177  for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1178 
1179  const double alpha =
1180  alpha_constant * t_w * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
1181 
1182  size_t rr = 0;
1183  for (; rr != OpBase::nbRows / SPACE_DIM; ++rr) {
1184  auto t_mat_diag = getFTensor1FromArrayDiag<SPACE_DIM, SPACE_DIM>(
1185  this->locMat, SPACE_DIM * rr);
1186  const double t_row_div_base = t_row_diff_base(i, i);
1187  auto t_col_base = col_data.getFTensor0N(gg, 0);
1188 
1189  for (size_t cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1190  t_mat_diag(i) += alpha * t_row_div_base * t_col_base;
1191  if constexpr (CoordSys == CYLINDRICAL) {
1192  t_mat_diag(i) += t_row_base(0) * (alpha / t_coords(0)) * t_col_base;
1193  }
1194  ++t_col_base;
1195  ++t_mat_diag;
1196  }
1197  ++t_row_base;
1198  ++t_row_diff_base;
1199  }
1200  for (; rr < nb_base_functions; ++rr) {
1201  ++t_row_diff_base;
1202  ++t_row_base;
1203  }
1204 
1205  ++t_w;
1206  ++t_coords;
1207  }
1208 
1209  MoFEMFunctionReturn(0);
1210 }
1211 
1212 template <int SPACE_DIM, typename OpBase, CoordinateTypes COORDINATE_SYSTEM>
1213 MoFEMErrorCode
1214 OpMixScalarTimesDivImpl<SPACE_DIM, GAUSS, OpBase, COORDINATE_SYSTEM>::iNtegrate(
1215  EntitiesFieldData::EntData &row_data,
1216  EntitiesFieldData::EntData &col_data) {
1217  MoFEMFunctionBegin;
1218 
1219 #ifndef NDEBUG
1220  if (OpBase::locMat.size2() % SPACE_DIM)
1221  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1222  "Number of rows in matrix should be multiple of space dimensions");
1223 #endif
1224 
1225  // When we move to C++17 add if constexpr()
1226  if constexpr (COORDINATE_SYSTEM == POLAR || COORDINATE_SYSTEM == SPHERICAL)
1227  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1228  "%s coordiante not implemented",
1229  CoordinateTypesNames[COORDINATE_SYSTEM]);
1230 
1231  auto t_w = this->getFTensor0IntegrationWeight();
1232  auto t_coords = this->getFTensor1CoordsAtGaussPts();
1233  size_t nb_base_functions_row = row_data.getN().size2();
1234  auto t_row_base = row_data.getFTensor0N();
1235  const double vol = this->getMeasure();
1236  for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1237 
1238  const double alpha =
1239  alphaConstant(t_coords(0), t_coords(1), t_coords(2)) * t_w * vol;
1240 
1241  size_t rr = 0;
1242  auto t_m = getFTensor1FromPtr<SPACE_DIM>(OpBase::locMat.data().data());
1243 
1244  // When we move to C++17 add if constexpr()
1245  if constexpr (COORDINATE_SYSTEM == CARTESIAN) {
1246  for (; rr != OpBase::nbRows; ++rr) {
1247  const double r_val = alpha * t_row_base;
1248  auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1249  for (size_t cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1250  t_m(i) += r_val * t_col_diff_base(i);
1251  ++t_col_diff_base;
1252  ++t_m;
1253  }
1254  ++t_row_base;
1255  }
1256  }
1257 
1258  // When we move to C++17 add if constexpr()
1259  if constexpr (COORDINATE_SYSTEM == CYLINDRICAL) {
1260  for (; rr != OpBase::nbRows; ++rr) {
1261  const double r_val = alpha * t_row_base;
1262  auto t_col_base = col_data.getFTensor0N(gg, 0);
1263  auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1264  for (size_t cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1265  t_m(i) += r_val * t_col_diff_base(i);
1266  t_m(0) += (r_val / t_coords(0)) * t_col_base;
1267  ++t_col_base;
1268  ++t_col_diff_base;
1269  ++t_m;
1270  }
1271  ++t_row_base;
1272  }
1273  }
1274 
1275  for (; rr < nb_base_functions_row; ++rr)
1276  ++t_row_base;
1277 
1278  ++t_w;
1279  ++t_coords;
1280  }
1281 
1282  MoFEMFunctionReturn(0);
1283 }
1284 
1285 template <int SPACE_DIM, typename OpBase>
1286 MoFEMErrorCode
1287 OpMixVectorTimesGradImpl<3, SPACE_DIM, SPACE_DIM, GAUSS, OpBase>::iNtegrate(
1288  EntitiesFieldData::EntData &row_data,
1289  EntitiesFieldData::EntData &col_data) {
1290  MoFEMFunctionBegin;
1291 
1292  auto t_w = this->getFTensor0IntegrationWeight();
1293 
1294  size_t nb_base_functions = row_data.getN().size2() / 3;
1295  auto t_row_base = row_data.getFTensor1N<3>();
1296  auto &mat = this->locMat;
1297  const double alpha_constant = alphaConstant();
1298  for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1299 
1300  const double alpha = alpha_constant * this->getMeasure() * t_w;
1301 
1302  size_t rr = 0;
1303  for (; rr != OpBase::nbRows; ++rr) {
1304  auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1305  for (size_t cc = 0; cc != OpBase::nbCols; ++cc) {
1306  mat(rr, cc) += alpha * t_row_base(i) * t_col_diff_base(i);
1307  ++t_col_diff_base;
1308  }
1309  ++t_row_base;
1310  }
1311  for (; rr < nb_base_functions; ++rr)
1312  ++t_row_base;
1313 
1314  ++t_w;
1315  }
1316 
1317  MoFEMFunctionReturn(0);
1318 }
1319 
1320 template <int SPACE_DIM, typename OpBase>
1321 MoFEMErrorCode
1322 OpMixVectorTimesGradImpl<1, SPACE_DIM, SPACE_DIM, GAUSS, OpBase>::iNtegrate(
1323  EntitiesFieldData::EntData &row_data,
1324  EntitiesFieldData::EntData &col_data) {
1325  MoFEMFunctionBegin;
1326 
1327  auto t_w = this->getFTensor0IntegrationWeight();
1328 
1329  size_t nb_base_functions = row_data.getN().size2();
1330  auto t_row_base = row_data.getFTensor0N();
1331 
1332  auto get_t_vec = [&](const int rr) {
1333  std::array<double *, SPACE_DIM> ptrs;
1334  for (auto i = 0; i != SPACE_DIM; ++i)
1335  ptrs[i] = &OpBase::locMat(rr + i, 0);
1336  return FTensor::Tensor1<FTensor::PackPtr<double *, 1>, SPACE_DIM>(ptrs);
1337  };
1338 
1339  const double alpha_constant = alphaConstant();
1340  for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1341 
1342  const double alpha = alpha_constant * this->getMeasure() * t_w;
1343 
1344  size_t rr = 0;
1345  for (; rr != OpBase::nbRows / SPACE_DIM; ++rr) {
1346  auto t_vec = get_t_vec(SPACE_DIM * rr);
1347  auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1348  for (size_t cc = 0; cc != OpBase::nbCols; ++cc) {
1349  t_vec(i) += alpha * t_row_base * t_col_diff_base(i);
1350  ++t_col_diff_base;
1351  ++t_vec;
1352  }
1353  ++t_row_base;
1354  }
1355  for (; rr < nb_base_functions; ++rr)
1356  ++t_row_base;
1357 
1358  ++t_w;
1359  }
1360 
1361  MoFEMFunctionReturn(0);
1362 }
1363 
1364 template <int SPACE_DIM, typename OpBase>
1365 MoFEMErrorCode OpMixTensorTimesGradImpl<SPACE_DIM, GAUSS, OpBase>::iNtegrate(
1366  EntitiesFieldData::EntData &row_data,
1367  EntitiesFieldData::EntData &col_data) {
1368  MoFEMFunctionBegin;
1369 
1370  auto t_w = this->getFTensor0IntegrationWeight();
1371  auto t_coords = this->getFTensor1CoordsAtGaussPts();
1372 
1373  size_t nb_base_functions = row_data.getN().size2() / 3;
1374  auto t_row_base = row_data.getFTensor1N<3>();
1375  const double alpha_constant = alphaConstant() * this->getMeasure();
1376  for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1377 
1378  const double alpha =
1379  alpha_constant * betaCoeff(t_coords(0), t_coords(1), t_coords(2)) * t_w;
1380 
1381  size_t rr = 0;
1382  for (; rr != OpBase::nbRows / SPACE_DIM; ++rr) {
1383  auto t_mat_diag = getFTensor1FromArrayDiag<SPACE_DIM, SPACE_DIM>(
1384  this->locMat, SPACE_DIM * rr);
1385  auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1386 
1387  for (size_t cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1388  t_mat_diag(i) += alpha * t_row_base(j) * t_col_diff_base(j);
1389  ++t_col_diff_base;
1390  ++t_mat_diag;
1391  }
1392 
1393  ++t_row_base;
1394  }
1395  for (; rr < nb_base_functions; ++rr)
1396  ++t_row_base;
1397 
1398  ++t_w;
1399  ++t_coords;
1400  }
1401 
1402  MoFEMFunctionReturn(0);
1403 }
1404 
1405 } // namespace MoFEM
1406 
1407 #endif //__BILINEAR_FORMS_INTEGRATORS_IMPL_HPP__
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::OpMixTensorTimesGradImpl< SPACE_DIM, GAUSS, OpBase >::j
FTensor::Index< 'j', SPACE_DIM > j
summit Index
Definition: BiLinearFormsIntegratorsImpl.hpp:414
MoFEM::OpMixTensorTimesGradImpl< SPACE_DIM, GAUSS, OpBase >::OpMixTensorTimesGradImpl
OpMixTensorTimesGradImpl(const std::string row_field_name, const std::string col_field_name, ConstantFun alpha_fun, ScalarFun beta_coeff, const bool assemble_transpose, const bool only_transpose=false)
Definition: BiLinearFormsIntegratorsImpl.hpp:401
MoFEM::OpBaseImpl::nbIntegrationPts
int nbIntegrationPts
number of integration points
Definition: FormsIntegrators.hpp:238
FTensor::Christof
Definition: Christof_value.hpp:9
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
MoFEM::OpMassCacheImpl< 1, FIELD_DIM, I, OpBase >::scalarBeta
double scalarBeta
Definition: BiLinearFormsIntegratorsImpl.hpp:177
MoFEM::EntitiesFieldData::EntData::getFieldEntities
const VectorFieldEntities & getFieldEntities() const
get field entities
Definition: EntitiesFieldData.hpp:1277
MoFEM::OpBaseImpl::locMat
MatrixDouble locMat
local entity block matrix
Definition: FormsIntegrators.hpp:249
MoFEM::OpMassCacheImpl< 1, FIELD_DIM, I, OpBase >::OpMassCacheImpl
OpMassCacheImpl(CacheMatsTypeType cache, const std::string row_field_name, const std::string col_field_name, const double beta, boost::shared_ptr< Range > ents_ptr=nullptr)
Definition: BiLinearFormsIntegratorsImpl.hpp:165
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
CARTESIAN
@ CARTESIAN
Definition: definitions.h:128
MoFEM::OpMassImpl< 1, 1, GAUSS, OpBase >::OpMassImpl
OpMassImpl(const std::string row_field_name, const std::string col_field_name, ScalarFun beta=[](double, double, double) constexpr { return 1;}, boost::shared_ptr< Range > ents_ptr=nullptr, boost::shared_ptr< MatrixDouble > cache_mat=nullptr)
Definition: BiLinearFormsIntegratorsImpl.hpp:69
MoFEM::OpGradSymTensorGradImpl< 1, SPACE_DIM, SPACE_DIM, S, GAUSS, OpBase >::OpGradSymTensorGradImpl
OpGradSymTensorGradImpl(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< MatrixDouble > mat_D, boost::shared_ptr< Range > ents_ptr=nullptr, ScalarFun beta=scalar_fun_one)
Definition: BiLinearFormsIntegratorsImpl.hpp:189
MoFEM::OpGradGradImpl< 1, 1, SPACE_DIM, GAUSS, OpBase >::OpGradGradImpl
OpGradGradImpl(const std::string row_field_name, const std::string col_field_name, ScalarFun beta=scalar_fun_one, boost::shared_ptr< Range > ents_ptr=nullptr)
Definition: BiLinearFormsIntegratorsImpl.hpp:28
MoFEM::OpMassImpl< 3, FIELD_DIM, GAUSS, OpBase >::OpMassImpl
OpMassImpl(const std::string row_field_name, const std::string col_field_name, ScalarFun beta=scalar_fun_one, boost::shared_ptr< Range > ents_ptr=nullptr)
Definition: BiLinearFormsIntegratorsImpl.hpp:110
MoFEM::OpGradGradSymTensorGradGradImpl< 1, 1, SPACE_DIM, S, GAUSS, OpBase >::i
FTensor::Index< 'i', SPACE_DIM > i
summit Index
Definition: BiLinearFormsIntegratorsImpl.hpp:236
MoFEM::OpMassImpl< 3, 9, GAUSS, OpBase >::OpMassImpl
OpMassImpl(const std::string row_field_name, const std::string col_field_name, ScalarFun beta=scalar_fun_one, boost::shared_ptr< Range > ents_ptr=nullptr)
Definition: BiLinearFormsIntegratorsImpl.hpp:144
MoFEM::OpMixDivTimesVecImpl< SPACE_DIM, GAUSS, OpBase, CoordSys >::OpMixDivTimesVecImpl
OpMixDivTimesVecImpl(const std::string row_field_name, const std::string col_field_name, ConstantFun alpha_fun, const bool assemble_transpose, const bool only_transpose=false)
Definition: BiLinearFormsIntegratorsImpl.hpp:284
MoFEM::OpMixDivTimesVecImpl< SPACE_DIM, GAUSS, OpBase, CoordSys >::i
FTensor::Index< 'i', SPACE_DIM > i
summit Index
Definition: BiLinearFormsIntegratorsImpl.hpp:305
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpMixDivTimesScalarImpl< SPACE_DIM, GAUSS, OpBase >::i
FTensor::Index< 'i', SPACE_DIM > i
summit Index
Definition: BiLinearFormsIntegratorsImpl.hpp:270
MoFEM::ConstantFun
boost::function< double()> ConstantFun
Constant function type.
Definition: FormsIntegrators.hpp:166
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
OpBase
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
BASE_DIM
constexpr int BASE_DIM
Definition: dg_projection.cpp:15
CoordinateTypes
CoordinateTypes
Coodinate system.
Definition: definitions.h:127
MoFEM::OpGradSymTensorGradImpl< 1, SPACE_DIM, SPACE_DIM, S, GAUSS, OpBase >::matD
boost::shared_ptr< MatrixDouble > matD
Definition: BiLinearFormsIntegratorsImpl.hpp:201
FIELD_DIM
constexpr int FIELD_DIM
Definition: child_and_parent.cpp:15
MoFEM::OpGradSymTensorGradImpl
Definition: BiLinearFormsIntegratorsImpl.hpp:183
I
constexpr IntegrationType I
Definition: operators_tests.cpp:31
MoFEM::scalar_fun_one
double scalar_fun_one(const double, const double, const double)
Definition: FormsIntegrators.hpp:146
MoFEM::OpMassCacheImpl< 1, FIELD_DIM, I, OpBase >::cacheLocMats
CacheMatsTypeType cacheLocMats
Definition: BiLinearFormsIntegratorsImpl.hpp:178
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1502
OpBase
MoFEM::OpMixScalarTimesDivImpl< SPACE_DIM, GAUSS, OpBase, COORDINATE_SYSTEM >::OpMixScalarTimesDivImpl
OpMixScalarTimesDivImpl(const std::string row_field_name, const std::string col_field_name, ScalarFun alpha_fun=[](double, double, double) constexpr { return 1;}, const bool assemble_transpose=false, const bool only_transpose=false)
Definition: BiLinearFormsIntegratorsImpl.hpp:321
MoFEM::OpGradSymTensorGradImpl< 1, SPACE_DIM, SPACE_DIM, S, GAUSS, OpBase >::i
FTensor::Index< 'i', SPACE_DIM > i
summit Index
Definition: BiLinearFormsIntegratorsImpl.hpp:188
MoFEM::OpGradGradImpl
Definition: BiLinearFormsIntegratorsImpl.hpp:23
MoFEM::EntitiesFieldData::EntData::getFTensor2DiffN
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
Definition: EntitiesFieldData.cpp:660
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:178
MoFEM::OpGradGradImpl< 1, 1, SPACE_DIM, GAUSS, OpBase >::i
FTensor::Index< 'i', SPACE_DIM > i
summit Index
Definition: BiLinearFormsIntegratorsImpl.hpp:27
CoordinateTypesNames
const static char *const CoordinateTypesNames[]
Coordinate system names.
Definition: definitions.h:121
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
FTensor::Tensor3
Definition: Tensor3_value.hpp:12
MoFEM::OpMixVectorTimesGradImpl< 3, SPACE_DIM, SPACE_DIM, GAUSS, OpBase >::alphaConstant
ConstantFun alphaConstant
Definition: BiLinearFormsIntegratorsImpl.hpp:359
MoFEM::OpGradGradSymTensorGradGradImpl< 1, 1, SPACE_DIM, S, GAUSS, OpBase >::OpGradGradSymTensorGradGradImpl
OpGradGradSymTensorGradGradImpl(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< MatrixDouble > mat_D, boost::shared_ptr< Range > ents_ptr=nullptr)
Definition: BiLinearFormsIntegratorsImpl.hpp:237
MoFEM::OpMassImpl< 1, FIELD_DIM, GAUSS, OpBase >::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: BiLinearFormsIntegratorsImpl.hpp:101
MoFEM::OpGradTensorGradImpl< 1, SPACE_DIM, SPACE_DIM, S, GAUSS, OpBase >::OpGradTensorGradImpl
OpGradTensorGradImpl(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< MatrixDouble > mat_D, ScalarFun beta=scalar_fun_one)
Definition: BiLinearFormsIntegratorsImpl.hpp:215
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::OpMassImpl< 3, 4, GAUSS, OpBase >::OpMassImpl
OpMassImpl(const std::string row_field_name, const std::string col_field_name, ScalarFun beta=scalar_fun_one, boost::shared_ptr< Range > ents_ptr=nullptr)
Definition: BiLinearFormsIntegratorsImpl.hpp:127
a
constexpr double a
Definition: approx_sphere.cpp:30
POLAR
@ POLAR
Definition: definitions.h:129
MoFEM::OpMixVectorTimesGradImpl< 3, SPACE_DIM, SPACE_DIM, GAUSS, OpBase >::i
FTensor::Index< 'i', SPACE_DIM > i
summit Index
Definition: BiLinearFormsIntegratorsImpl.hpp:358
MoFEM::OpMixTensorTimesGradImpl< SPACE_DIM, GAUSS, OpBase >::OpMixTensorTimesGradImpl
OpMixTensorTimesGradImpl(const std::string row_field_name, const std::string col_field_name, ConstantFun alpha_fun, const bool assemble_transpose, const bool only_transpose=false)
Definition: BiLinearFormsIntegratorsImpl.hpp:391
MoFEM::OpGradTensorGradImpl
Definition: BiLinearFormsIntegratorsImpl.hpp:209
double
MoFEM::OpMassImpl< 3, FIELD_DIM, GAUSS, OpBase >::betaCoeff
ScalarFun betaCoeff
Definition: BiLinearFormsIntegratorsImpl.hpp:120
MoFEM::OpMassImpl< 3, 9, GAUSS, OpBase >::betaCoeff
ScalarFun betaCoeff
Definition: BiLinearFormsIntegratorsImpl.hpp:154
MoFEM::OpMixDivTimesVecImpl
Definition: BiLinearFormsIntegratorsImpl.hpp:278
MoFEM::OpGradSymTensorGradImpl< 1, SPACE_DIM, SPACE_DIM, S, GAUSS, OpBase >::betaCoeff
ScalarFun betaCoeff
Definition: BiLinearFormsIntegratorsImpl.hpp:202
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1214
MoFEM::OpGradTensorGradImpl< 1, SPACE_DIM, SPACE_DIM, S, GAUSS, OpBase >::betaCoeff
ScalarFun betaCoeff
Definition: BiLinearFormsIntegratorsImpl.hpp:224
MoFEM::OpMixVectorTimesGradImpl< 3, SPACE_DIM, SPACE_DIM, GAUSS, OpBase >::OpMixVectorTimesGradImpl
OpMixVectorTimesGradImpl(const std::string row_field_name, const std::string col_field_name, ConstantFun alpha_fun, const bool assemble_transpose=false, const bool only_transpose=false)
Definition: BiLinearFormsIntegratorsImpl.hpp:346
MoFEM::OpGradGradImpl< 1, FIELD_DIM, SPACE_DIM, GAUSS, OpBase >::i
FTensor::Index< 'i', SPACE_DIM > i
summit Index
Definition: BiLinearFormsIntegratorsImpl.hpp:46
SPHERICAL
@ SPHERICAL
Definition: definitions.h:131
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
MoFEM::OpBaseImpl::colSide
int colSide
column side number
Definition: FormsIntegrators.hpp:242
MoFEM::OpGradGradImpl< 1, FIELD_DIM, SPACE_DIM, GAUSS, OpBase >::OpGradGradImpl
OpGradGradImpl(const std::string row_field_name, const std::string col_field_name, ScalarFun beta=scalar_fun_one, boost::shared_ptr< Range > ents_ptr=nullptr)
Definition: BiLinearFormsIntegratorsImpl.hpp:47
MoFEM::OpGradGradImpl< 1, 1, SPACE_DIM, GAUSS, OpBase >::betaCoeff
ScalarFun betaCoeff
Definition: BiLinearFormsIntegratorsImpl.hpp:39
MoFEM::ScalarFun
boost::function< double(const double, const double, const double)> ScalarFun
Scalar function type.
Definition: FormsIntegrators.hpp:144
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::OpGradTensorGradImpl< 1, SPACE_DIM, SPACE_DIM, S, GAUSS, OpBase >::i
FTensor::Index< 'i', SPACE_DIM > i
summit Index
Definition: BiLinearFormsIntegratorsImpl.hpp:214
MoFEM::OpMixScalarTimesDivImpl
Definition: BiLinearFormsIntegratorsImpl.hpp:316
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
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::rowSide
int rowSide
row side number
Definition: FormsIntegrators.hpp:241
MoFEM::OpBaseImpl::nbCols
int nbCols
number if dof on column
Definition: FormsIntegrators.hpp:237
MoFEM::OpMixVectorTimesGradImpl
Definition: BiLinearFormsIntegratorsImpl.hpp:341
MoFEM::OpMassImpl
Definition: BiLinearFormsIntegratorsImpl.hpp:64
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
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
CYLINDRICAL
@ CYLINDRICAL
Definition: definitions.h:130
MoFEM::OpMixDivTimesScalarImpl< SPACE_DIM, GAUSS, OpBase >::alphaConstant
ConstantFun alphaConstant
Definition: BiLinearFormsIntegratorsImpl.hpp:271
MoFEM::OpMixTensorTimesGradImpl
Definition: BiLinearFormsIntegratorsImpl.hpp:386
MoFEM::OpMassImpl< 1, 1, GAUSS, OpBase >::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: BiLinearFormsIntegratorsImpl.hpp:85
MoFEM::OpMixVectorTimesGradImpl< 1, SPACE_DIM, SPACE_DIM, GAUSS, OpBase >::OpMixVectorTimesGradImpl
OpMixVectorTimesGradImpl(const std::string row_field_name, const std::string col_field_name, ConstantFun alpha_fun, const bool assemble_transpose=false, const bool only_transpose=false)
Definition: BiLinearFormsIntegratorsImpl.hpp:367
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: EntitiesFieldData.hpp:1318
MoFEM::OpMixVectorTimesGradImpl< 1, SPACE_DIM, SPACE_DIM, GAUSS, OpBase >::alphaConstant
ConstantFun alphaConstant
Definition: BiLinearFormsIntegratorsImpl.hpp:380
MoFEM::OpMassCacheImpl
Definition: BiLinearFormsIntegratorsImpl.hpp:160
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::EntitiesFieldData::EntData::getFTensor1N
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
Definition: EntitiesFieldData.cpp:640
MoFEM::OpMixVectorTimesGradImpl< 1, SPACE_DIM, SPACE_DIM, GAUSS, OpBase >::i
FTensor::Index< 'i', SPACE_DIM > i
summit Index
Definition: BiLinearFormsIntegratorsImpl.hpp:379
MoFEM::OpMixDivTimesScalarImpl
Definition: BiLinearFormsIntegratorsImpl.hpp:254
MoFEM::OpGradGradSymTensorGradGradImpl
Definition: BiLinearFormsIntegratorsImpl.hpp:231
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::OpGradGradImpl< 1, FIELD_DIM, SPACE_DIM, GAUSS, OpBase >::betaCoeff
ScalarFun betaCoeff
Definition: BiLinearFormsIntegratorsImpl.hpp:58
MoFEM::OpMixDivTimesScalarImpl< SPACE_DIM, GAUSS, OpBase >::OpMixDivTimesScalarImpl
OpMixDivTimesScalarImpl(const std::string row_field_name, const std::string col_field_name, ConstantFun alpha_fun, const bool assemble_transpose=false, const bool only_transpose=false)
Definition: BiLinearFormsIntegratorsImpl.hpp:258
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
MoFEM::OpGradGradSymTensorGradGradImpl< 1, 1, SPACE_DIM, S, GAUSS, OpBase >::matD
boost::shared_ptr< MatrixDouble > matD
Definition: BiLinearFormsIntegratorsImpl.hpp:248
MoFEM::OpMixScalarTimesDivImpl< SPACE_DIM, GAUSS, OpBase, COORDINATE_SYSTEM >::alphaConstant
ScalarFun alphaConstant
Definition: BiLinearFormsIntegratorsImpl.hpp:334
MoFEM::CacheMatsTypeType
std::map< std::pair< int, int >, boost::shared_ptr< MatrixDouble > > CacheMatsTypeType
Definition: BiLinearFormsIntegratorsImpl.hpp:19
MoFEM::OpBaseImpl::rowType
EntityType rowType
row type
Definition: FormsIntegrators.hpp:243
MoFEM::OpMixTensorTimesGradImpl< SPACE_DIM, GAUSS, OpBase >::i
FTensor::Index< 'i', SPACE_DIM > i
summit Index
Definition: BiLinearFormsIntegratorsImpl.hpp:413
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::OpMixDivTimesVecImpl< SPACE_DIM, GAUSS, OpBase, CoordSys >::OpMixDivTimesVecImpl
OpMixDivTimesVecImpl(const std::string row_field_name, const std::string col_field_name, ConstantFun alpha_fun, ScalarFun beta_fun, const bool assemble_transpose, const bool only_transpose=false)
Definition: BiLinearFormsIntegratorsImpl.hpp:294
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEM::OpGradTensorGradImpl< 1, SPACE_DIM, SPACE_DIM, S, GAUSS, OpBase >::matD
boost::shared_ptr< MatrixDouble > matD
Definition: BiLinearFormsIntegratorsImpl.hpp:223
MoFEM::OpMassImpl< 3, 4, GAUSS, OpBase >::betaCoeff
ScalarFun betaCoeff
Definition: BiLinearFormsIntegratorsImpl.hpp:137
MoFEM::OpMixScalarTimesDivImpl< SPACE_DIM, GAUSS, OpBase, COORDINATE_SYSTEM >::i
FTensor::Index< 'i', SPACE_DIM > i
summit Index
Definition: BiLinearFormsIntegratorsImpl.hpp:333
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MoFEM::OpBaseImpl::colType
EntityType colType
column type
Definition: FormsIntegrators.hpp:244
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
Generated by Doxygen 1.8.17 and hosted at University of Glasgow