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