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