v0.13.1
BiLinearFormsIntegrators.hpp
Go to the documentation of this file.
1/** \file BiLinearFormsIntegrators.hpp
2 * \brief Bilinear forms integrators
3 * \ingroup mofem_form
4
5 \todo SSome 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
14
15#ifndef __BILINEAR_FORMS_INTEGRATORS_HPP__
16#define __BILINEAR_FORMS_INTEGRATORS_HPP__
17
18namespace MoFEM {
19
20template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, IntegrationType I,
21 typename OpBase>
23
24template <int SPACE_DIM, typename OpBase>
25struct OpGradGradImpl<1, 1, SPACE_DIM, GAUSS, OpBase> : public OpBase {
27 OpGradGradImpl(const std::string row_field_name,
28 const std::string col_field_name, ScalarFun beta,
29 boost::shared_ptr<Range> ents_ptr = nullptr)
30 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
31 betaCoeff(beta) {
32 if (row_field_name == col_field_name)
33 this->sYmm = true;
34 }
35
36protected:
38 boost::shared_ptr<Range> entsPtr;
41};
42
43template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
46 OpGradGradImpl(const std::string row_field_name,
47 const std::string col_field_name, ScalarFun beta,
48 boost::shared_ptr<Range> ents_ptr = nullptr)
49 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
50 betaCoeff(beta) {
51 if (row_field_name == col_field_name)
52 this->sYmm = true;
53 }
54
55protected:
57 boost::shared_ptr<Range> entsPtr;
60};
61
62template <int BASE_DIM, int FIELD_DIM, IntegrationType I, typename OpBase>
63struct OpMassImpl {};
64
65template <typename OpBase>
66struct OpMassImpl<1, 1, GAUSS, OpBase> : public OpBase {
67
68 OpMassImpl(const std::string row_field_name, const std::string col_field_name,
69 ScalarFun beta, boost::shared_ptr<Range> ents_ptr = nullptr)
70 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
71 betaCoeff(beta), entsPtr(ents_ptr) {
72 if (row_field_name == col_field_name)
73 this->sYmm = true;
74 }
75
76protected:
78 boost::shared_ptr<Range> entsPtr;
81};
82
83template <int FIELD_DIM, typename OpBase>
85 : public OpMassImpl<1, 1, GAUSS, OpBase> {
86 using OpMassImpl<1, 1, GAUSS, OpBase>::OpMassImpl;
87
88protected:
89 using OpMassImpl<1, 1, GAUSS, OpBase>::betaCoeff;
90 using OpMassImpl<1, 1, GAUSS, OpBase>::entsPtr;
93};
94
95template <int BASE_DIM, typename OpBase>
97
98 OpMassImpl(const std::string row_field_name, const std::string col_field_name,
99 ScalarFun beta, boost::shared_ptr<Range> ents_ptr = nullptr)
100 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
101 betaCoeff(beta), entsPtr(ents_ptr) {
102 if (row_field_name == col_field_name)
103 this->sYmm = true;
104 }
105
106protected:
108 boost::shared_ptr<Range> entsPtr;
109 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
111};
112
113template <typename OpBase>
114struct OpMassImpl<3, 9, GAUSS, OpBase> : public OpBase {
115 OpMassImpl(const std::string row_field_name, const std::string col_field_name,
116 ScalarFun beta, boost::shared_ptr<Range> ents_ptr = nullptr)
117 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
118 betaCoeff(beta), entsPtr(ents_ptr) {
119 if (row_field_name == col_field_name)
120 this->sYmm = true;
121 }
122
123protected:
125 boost::shared_ptr<Range> entsPtr;
126 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
128};
129
130template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S, IntegrationType I,
131 typename OpBase>
133
134template <int SPACE_DIM, int S, typename OpBase>
136 : public OpBase {
138 OpGradSymTensorGradImpl(const std::string row_field_name,
139 const std::string col_field_name,
140 boost::shared_ptr<MatrixDouble> mat_D,
141 boost::shared_ptr<Range> ents_ptr = nullptr)
142 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL), matD(mat_D),
143 entsPtr(ents_ptr) {
144 if (row_field_name == col_field_name)
145 this->sYmm = true;
146 }
147
148protected:
149 boost::shared_ptr<MatrixDouble> matD;
150 boost::shared_ptr<Range> entsPtr;
151 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
153};
154
155template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S, IntegrationType I,
156 typename OpBase>
158
159template <int SPACE_DIM, int S, typename OpBase>
161 : public OpBase {
163 OpGradTensorGradImpl(const std::string row_field_name,
164 const std::string col_field_name,
165 boost::shared_ptr<MatrixDouble> mat_D)
166 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL), matD(mat_D) {}
167
168protected:
169 boost::shared_ptr<MatrixDouble> matD;
170 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
172};
173
174template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S, IntegrationType I,
175 typename OpBase>
177
178template <int SPACE_DIM, int S, typename OpBase>
180 : public OpBase {
182 OpGradGradSymTensorGradGradImpl(const std::string row_field_name,
183 const std::string col_field_name,
184 boost::shared_ptr<MatrixDouble> mat_D,
185 boost::shared_ptr<Range> ents_ptr = nullptr)
186 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL), matD(mat_D) {
187 if (row_field_name == col_field_name)
188 this->sYmm = true;
189 }
190
191protected:
192 boost::shared_ptr<MatrixDouble> matD;
193 boost::shared_ptr<Range> entsPtr;
194 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
196};
197
198template <int SPACE_DIM, IntegrationType I, typename OpBase>
200
201template <int SPACE_DIM, typename OpBase>
203 OpMixDivTimesScalarImpl(const std::string row_field_name,
204 const std::string col_field_name,
205 ConstantFun alpha_fun,
206 const bool assemble_transpose = false,
207 const bool only_transpose = false)
208 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
209 alphaConstant(alpha_fun) {
210 this->assembleTranspose = assemble_transpose;
211 this->onlyTranspose = only_transpose;
212 }
213
214protected:
217 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
219};
220
221template <int SPACE_DIM, IntegrationType I, typename OpBase>
223
224template <int SPACE_DIM, typename OpBase>
226 OpMixDivTimesVecImpl(const std::string row_field_name,
227 const std::string col_field_name, ConstantFun alpha_fun,
228 const bool assemble_transpose = false,
229 const bool only_transpose = false)
230 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
231 alphaConstant(alpha_fun) {
232 this->assembleTranspose = assemble_transpose;
233 this->onlyTranspose = only_transpose;
234 }
235
236protected:
239 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
241};
242
243template <int SPACE_DIM, IntegrationType I, typename OpBase,
244 CoordinateTypes COORDINATE_SYSTEM>
246
247template <int SPACE_DIM, typename OpBase, CoordinateTypes COORDINATE_SYSTEM>
248struct OpMixScalarTimesDivImpl<SPACE_DIM, GAUSS, OpBase, COORDINATE_SYSTEM>
249 : public OpBase {
250 OpMixScalarTimesDivImpl(const std::string row_field_name,
251 const std::string col_field_name, ScalarFun alpha_fun,
252 const bool assemble_transpose = false,
253 const bool only_transpose = false)
254 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
255 alphaConstant(alpha_fun) {
256 this->assembleTranspose = assemble_transpose;
257 this->onlyTranspose = only_transpose;
258 this->sYmm = false;
259 }
260
261protected:
264 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
266};
267
268template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, IntegrationType I,
269 typename OpBase>
271
272template <int SPACE_DIM, typename OpBase>
274 : public OpBase {
275 OpMixVectorTimesGradImpl(const std::string row_field_name,
276 const std::string col_field_name,
277 ConstantFun alpha_fun,
278 const bool assemble_transpose = false,
279 const bool only_transpose = false)
280 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
281 alphaConstant(alpha_fun) {
282 this->assembleTranspose = assemble_transpose;
283 this->onlyTranspose = only_transpose;
284 }
285
286protected:
289 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
291};
292
293template <int SPACE_DIM, typename OpBase>
295 : public OpBase {
296 OpMixVectorTimesGradImpl(const std::string row_field_name,
297 const std::string col_field_name,
298 ConstantFun alpha_fun,
299 const bool assemble_transpose = false,
300 const bool only_transpose = false)
301 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
302 alphaConstant(alpha_fun) {
303 this->assembleTranspose = assemble_transpose;
304 this->onlyTranspose = only_transpose;
305 }
306
307protected:
310 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
312};
313
314template <int SPACE_DIM, IntegrationType I, typename OpBase>
316
317template <int SPACE_DIM, typename OpBase>
319 OpMixTensorTimesGradImpl(const std::string row_field_name,
320 const std::string col_field_name,
321 ConstantFun alpha_fun,
322 const bool assemble_transpose = false,
323 const bool only_transpose = false)
324 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
325 alphaConstant(alpha_fun) {
326 this->assembleTranspose = assemble_transpose;
327 this->onlyTranspose = only_transpose;
328 }
329
330protected:
334 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
336};
337
338template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, IntegrationType I,
339 typename OpBase>
341
342template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, IntegrationType I,
343 typename OpBase>
345
346template <int SPACE_DIM, typename OpBase>
348 : public OpBase {
350 const std::string field_name_row, const std::string field_name_col,
351 boost::shared_ptr<MatrixDouble> y_grad_ptr,
352 ConstantFun alpha_fun = []() { return 1; })
353 : OpBase(field_name_row, field_name_col, OpBase::OPROWCOL),
354 yGradPtr(y_grad_ptr), alphaConstant(alpha_fun) {
355
356 this->assembleTranspose = false;
357 this->onlyTranspose = false;
358 this->sYmm = false;
359 }
360
361protected:
362 boost::shared_ptr<MatrixDouble> yGradPtr;
364 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
366};
367
368template <int SPACE_DIM, typename OpBase>
370 : public OpBase {
372 const std::string field_name_row, const std::string field_name_col,
373 boost::shared_ptr<MatrixDouble> u_ptr,
374 ConstantFun alpha_fun = []() { return 1; })
375 : OpBase(field_name_row, field_name_col, OpBase::OPROWCOL), uPtr(u_ptr),
376 alphaConstant(alpha_fun) {
377
378 this->assembleTranspose = false;
379 this->onlyTranspose = false;
380 this->sYmm = false;
381 }
382
383protected:
385 boost::shared_ptr<MatrixDouble> uPtr;
386 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
388};
389
390template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
392 : public OpBase {
394 const std::string field_name_row, const std::string field_name_col,
395 boost::shared_ptr<MatrixDouble> y_grad_ptr,
396 ConstantFun alpha_fun = []() { return 1; })
397 : OpBase(field_name_row, field_name_col, OpBase::OPROWCOL),
398 yGradPtr(y_grad_ptr), alphaConstant(alpha_fun) {
399
400 this->assembleTranspose = false;
401 this->onlyTranspose = false;
402 this->sYmm = false;
403 }
404
405protected:
407 boost::shared_ptr<MatrixDouble> yGradPtr;
408 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
410};
411
412template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
414 : public OpBase {
416 const std::string field_name_row, const std::string field_name_col,
417 boost::shared_ptr<MatrixDouble> u_ptr,
418 ConstantFun alpha_fun = []() { return 1; })
419 : OpBase(field_name_row, field_name_col, OpBase::OPROWCOL), uPtr(u_ptr),
420 alphaConstant(alpha_fun) {
421
422 this->assembleTranspose = false;
423 this->onlyTranspose = false;
424 this->sYmm = false;
425 }
426
427protected:
429 boost::shared_ptr<MatrixDouble> uPtr;
430 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
432};
433
434/**
435 * @brief Bilinear integrator form
436 * @ingroup mofem_forms
437 *
438 * @tparam EleOp
439 * @tparam A
440 * @tparam I
441 */
442template <typename EleOp>
443template <AssemblyType A>
444template <IntegrationType I>
446
447 /**
448 * @brief Integrate \f$(v_{,i},\beta(\mathbf{x}) u_{,j}))_\Omega\f$
449 * @ingroup mofem_forms
450 *
451 * @tparam SPACE_DIM
452 */
453 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM>
455 : public OpGradGradImpl<BASE_DIM, FIELD_DIM, SPACE_DIM, I, OpBase> {
457 OpBase>::OpGradGradImpl;
458 };
459
460 /**
461 * @brief Integrate \f$(v_i,\beta(\mathbf{x}) u_j)_\Omega\f$
462 * @ingroup mofem_forms
463 *
464 * @tparam
465 */
466 template <int BASE_DIM, int FIELD_DIM>
467 struct OpMass : public OpMassImpl<BASE_DIM, FIELD_DIM, I, OpBase> {
468 using OpMassImpl<BASE_DIM, FIELD_DIM, I, OpBase>::OpMassImpl;
469 };
470
471 /**
472 * @brief Integrate \f$(v_k,D_{ijkl} u_{,l})_\Omega\f$
473 *
474 * \note \f$D_{ijkl}\f$ is * tensor with minor & major symmetry
475 *
476 * @ingroup mofem_forms
477 *
478 * @tparam SPACE_DIM
479 * @tparam S
480 */
481 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S = 0>
482 struct OpGradSymTensorGrad
483 : public OpGradSymTensorGradImpl<BASE_DIM, FIELD_DIM, SPACE_DIM, S, I,
484 OpBase> {
486 OpBase>::OpGradSymTensorGradImpl;
487 };
488
489 /**
490 * @brief Integrate \f$(v_{,ij},D_{ijkl} u_{,kl})_\Omega\f$
491 *
492 * \note \f$D_{ijkl}\f$ is * tensor with no symmetries
493 *
494 * @ingroup mofem_forms
495 *
496 * @tparam SPACE_DIM
497 * @tparam S
498 */
499 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S = 0>
500 struct OpGradGradSymTensorGradGrad
501 : public OpGradGradSymTensorGradGradImpl<BASE_DIM, FIELD_DIM, SPACE_DIM,
502 S, I, OpBase> {
505 OpBase>::OpGradGradSymTensorGradGradImpl;
506 };
507
508 /**
509 * @brief Integrate \f$(v_{,j},D_{ijkl} u_{,l})_\Omega\f$
510 *
511 * \note \f$D_{ijkl}\f$ is * tensor with no symmetries
512 *
513 * @ingroup mofem_forms
514 *
515 * @tparam SPACE_DIM
516 * @tparam S
517 */
518 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S = 0>
519 struct OpGradTensorGrad
520 : public OpGradTensorGradImpl<BASE_DIM, FIELD_DIM, SPACE_DIM, S, I,
521 OpBase> {
523 OpBase>::OpGradTensorGradImpl;
524 };
525
526 /**
527 * @brief Integrate \f$(\lambda_{i,i},u)_\Omega\f$
528 *
529 * @tparam SPACE_DIM
530 */
531 template <int SPACE_DIM>
532 struct OpMixDivTimesScalar
533 : public OpMixDivTimesScalarImpl<SPACE_DIM, I, OpBase> {
535 OpBase>::OpMixDivTimesScalarImpl;
536 };
537
538 /**
539 * @brief Integrate \f$(\lambda_{ij,j},u_{i})_\Omega\f$
540 *
541 * @tparam SPACE_DIM
542 */
543 template <int SPACE_DIM>
544 struct OpMixDivTimesVec : public OpMixDivTimesVecImpl<SPACE_DIM, I, OpBase> {
545 using OpMixDivTimesVecImpl<SPACE_DIM, I, OpBase>::OpMixDivTimesVecImpl;
546 };
547
548 /**
549 * @brief Integrate \f$(\lambda,u_{i,i})_\Omega\f$
550 *
551 * @tparam SPACE_DIM
552 */
553 template <int SPACE_DIM, CoordinateTypes COORDINATE_SYSTEM = CARTESIAN>
555 : public OpMixScalarTimesDivImpl<SPACE_DIM, I, OpBase,
556 COORDINATE_SYSTEM> {
558 COORDINATE_SYSTEM>::OpMixScalarTimesDivImpl;
559 };
560
561 /**
562 * @brief Integrate \f$(\lambda_{i},u_{,j})_\Omega\f$
563 *
564 * @tparam SPACE_DIM
565 */
566 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM>
567 struct OpMixVectorTimesGrad
568 : public OpMixVectorTimesGradImpl<BASE_DIM, FIELD_DIM, SPACE_DIM, I,
569 OpBase> {
571 OpBase>::OpMixVectorTimesGradImpl;
572 };
573
574 /**
575 * @brief Integrate \f$(\lambda_{ij},u_{i,j})_\Omega\f$
576 *
577 * @tparam SPACE_DIM
578 */
579 template <int SPACE_DIM>
580 struct OpMixTensorTimesGrad
581 : public OpMixTensorTimesGradImpl<SPACE_DIM, I, OpBase> {
583 OpBase>::OpMixTensorTimesGradImpl;
584 };
585
586 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM>
587 struct OpConvectiveTermLhsDu
588 : public OpConvectiveTermLhsDuImpl<BASE_DIM, FIELD_DIM, SPACE_DIM, I,
589 OpBase> {
591 OpBase>::OpConvectiveTermLhsDuImpl;
592 };
593
594 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM>
595 struct OpConvectiveTermLhsDy
596 : public OpConvectiveTermLhsDyImpl<BASE_DIM, FIELD_DIM, SPACE_DIM, I,
597 OpBase> {
599 OpBase>::OpConvectiveTermLhsDyImpl;
600 };
601};
602
603template <int SPACE_DIM, typename OpBase>
606 EntitiesFieldData::EntData &col_data) {
608
609 if (entsPtr) {
610 if (entsPtr->find(OpBase::getFEEntityHandle()) == entsPtr->end())
612 }
613
614 // get element volume
615 const double vol = OpBase::getMeasure();
616 // get integration weights
617 auto t_w = OpBase::getFTensor0IntegrationWeight();
618 // get base function gradient on rows
619 auto t_row_grad = row_data.getFTensor1DiffN<SPACE_DIM>();
620 // get coordinate at integration points
621 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
622
623 // loop over integration points
624 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
625 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
626 // take into account Jacobian
627 const double alpha = t_w * beta;
628 // loop over ros base functions
629 int rr = 0;
630 for (; rr != OpBase::nbRows; rr++) {
631 // get column base functions gradient at gauss point gg
632 auto t_col_grad = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
633 // loop over columns
634 for (int cc = 0; cc != OpBase::nbCols; cc++) {
635 // calculate element of local matrix
636 OpBase::locMat(rr, cc) += alpha * (t_row_grad(i) * t_col_grad(i));
637 ++t_col_grad; // move to another gradient of base function
638 // on column
639 }
640 ++t_row_grad; // move to another element of gradient of base
641 // function on row
642 }
643 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
644 ++t_row_grad;
645
646 ++t_coords;
647 ++t_w; // move to another integration weight
648 }
650}
651
652template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
656 EntitiesFieldData::EntData &col_data) {
658
659 if (entsPtr) {
660 if (entsPtr->find(OpBase::getFEEntityHandle()) == entsPtr->end())
662 }
663
664 // get element volume
665 const double vol = OpBase::getMeasure();
666 // get integration weights
667 auto t_w = OpBase::getFTensor0IntegrationWeight();
668 // get base function gradient on rows
669 auto t_row_grad = row_data.getFTensor1DiffN<SPACE_DIM>();
670 // get coordinate at integration points
671 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
672
675
676 auto get_t_vec = [&](const int rr) {
677 std::array<double *, FIELD_DIM> ptrs;
678 for (auto i = 0; i != FIELD_DIM; ++i)
679 ptrs[i] = &OpBase::locMat(rr + i, i);
681 ptrs);
682 };
683
684 // loop over integration points
685 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
686 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
687 // take into account Jacobian
688 const double alpha = t_w * beta;
689 // loop over ros base functions
690 int rr = 0;
691 for (; rr != OpBase::nbRows / FIELD_DIM; rr++) {
692 // get diag vec
693 auto t_vec = get_t_vec(rr * FIELD_DIM);
694 // get column base functions gradient at gauss point gg
695 auto t_col_grad = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
696 // loop over columns
697 for (int cc = 0; cc != OpBase::nbCols / FIELD_DIM; cc++) {
698 // calculate element of local matrix
699 t_vec(j) += alpha * (t_row_grad(i) * t_col_grad(i));
700 ++t_col_grad; // move to another gradient of base function
701 // on column
702 ++t_vec;
703 }
704 ++t_row_grad; // move to another element of gradient of base
705 // function on row
706 }
707 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
708 ++t_row_grad;
709
710 ++t_coords;
711 ++t_w; // move to another integration weight
712 }
714}
715
716template <typename OpBase>
719 EntitiesFieldData::EntData &col_data) {
721 if (entsPtr) {
722 if (entsPtr->find(OpBase::getFEEntityHandle()) == entsPtr->end())
724 }
725
726#ifndef NDEBUG
727 auto log_error = [&]() {
728 MOFEM_LOG("SELF", Sev::error) << "Row side " << OpBase::rowSide << " "
729 << CN::EntityTypeName(OpBase::rowType);
730 MOFEM_LOG("SELF", Sev::error) << "Col side " << OpBase::colSide << " "
731 << CN::EntityTypeName(OpBase::colType);
732 };
733
734 if (row_data.getN().size2() < OpBase::nbRows) {
735 log_error();
736 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
737 "Wrong number of base functions on rows %d < %d",
738 row_data.getN().size2(), OpBase::nbRows);
739 }
740 if (col_data.getN().size2() < OpBase::nbCols) {
741 log_error();
742 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
743 "Wrong number of base functions on cols %d < %d",
744 col_data.getN().size2(), OpBase::nbCols);
745 }
746 if (row_data.getN().size1() != OpBase::nbIntegrationPts) {
747 log_error();
748 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
749 "Wrong number of integration points on rows %d != %d",
750 row_data.getN().size1(), OpBase::nbIntegrationPts);
751 }
752 if (col_data.getN().size1() != OpBase::nbIntegrationPts) {
753 log_error();
754 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
755 "Wrong number of integration points on cols %d != %d",
756 col_data.getN().size1(), OpBase::nbIntegrationPts);
757 }
758#endif
759
760 // get element volume
761 const double vol = OpBase::getMeasure();
762 // get integration weights
763 auto t_w = OpBase::getFTensor0IntegrationWeight();
764 // get base function gradient on rows
765 auto t_row_base = row_data.getFTensor0N();
766 // get coordinate at integration points
767 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
768 // loop over integration points
769 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
770 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
771 // take into account Jacobian
772 const double alpha = t_w * beta;
773 // loop over rows base functions
774 auto a_mat_ptr = &*OpBase::locMat.data().begin();
775 int rr = 0;
776 for (; rr != OpBase::nbRows; rr++) {
777 // get column base functions gradient at gauss point gg
778 auto t_col_base = col_data.getFTensor0N(gg, 0);
779 // loop over columns
780 for (int cc = 0; cc != OpBase::nbCols; cc++) {
781 // calculate element of local matrix
782 *a_mat_ptr += alpha * (t_row_base * t_col_base);
783 ++t_col_base;
784 ++a_mat_ptr;
785 }
786 ++t_row_base;
787 }
788 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
789 ++t_row_base;
790 ++t_coords;
791 ++t_w; // move to another integration weight
792 }
794};
795
796template <int FIELD_DIM, typename OpBase>
799 EntitiesFieldData::EntData &col_data) {
801 if (entsPtr) {
802 if (entsPtr->find(OpBase::getFEEntityHandle()) == entsPtr->end())
804 }
805 // get element volume
806 const double vol = OpBase::getMeasure();
807 // get integration weights
808 auto t_w = OpBase::getFTensor0IntegrationWeight();
809 // get base function gradient on rows
810 auto t_row_base = row_data.getFTensor0N();
811 // get coordinate at integration points
812 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
813
815 auto get_t_vec = [&](const int rr) {
816 std::array<double *, FIELD_DIM> ptrs;
817 for (auto i = 0; i != FIELD_DIM; ++i)
818 ptrs[i] = &OpBase::locMat(rr + i, i);
820 ptrs);
821 };
822
823 // loop over integration points
824 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
825 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
826 // take into account Jacobian
827 const double alpha = t_w * beta;
828 // loop over rows base functions
829 int rr = 0;
830 for (; rr != OpBase::nbRows / FIELD_DIM; rr++) {
831 // get column base functions gradient at gauss point gg
832 auto t_col_base = col_data.getFTensor0N(gg, 0);
833 // get mat vec
834 auto t_vec = get_t_vec(FIELD_DIM * rr);
835 // loop over columns
836 for (int cc = 0; cc != OpBase::nbCols / FIELD_DIM; cc++) {
837 // calculate element of local matrix
838 t_vec(i) += alpha * (t_row_base * t_col_base);
839 ++t_col_base;
840 ++t_vec;
841 }
842 ++t_row_base;
843 }
844 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
845 ++t_row_base;
846 ++t_coords;
847 ++t_w; // move to another integration weight
848 }
850};
851
852template <int BASE_DIM, typename OpBase>
855 EntitiesFieldData::EntData &col_data) {
858 if (entsPtr) {
859 if (entsPtr->find(OpBase::getFEEntityHandle()) == entsPtr->end())
861 }
862 size_t nb_base_functions = row_data.getN().size2() / BASE_DIM;
863 // // get element volume
864 const double vol = OpBase::getMeasure();
865 // get integration weights
866 auto t_w = OpBase::getFTensor0IntegrationWeight();
867 // get base function gradient on rows
868 auto t_row_base = row_data.getFTensor1N<BASE_DIM>();
869 // get coordinate at integration points
870 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
871 // loop over integration points
872 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
873 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
874 // take into account Jacobian
875 const double alpha = t_w * beta;
876 // loop over rows base functions
877 auto a_mat_ptr = &*OpBase::locMat.data().begin();
878 int rr = 0;
879 for (; rr != OpBase::nbRows; rr++) {
880 // get column base functions gradient at gauss point gg
881 auto t_col_base = col_data.getFTensor1N<BASE_DIM>(gg, 0);
882 // loop over columns
883 for (int cc = 0; cc != OpBase::nbCols; cc++) {
884 // calculate element of local matrix
885 (*a_mat_ptr) += alpha * (t_row_base(i) * t_col_base(i));
886 ++t_col_base;
887 ++a_mat_ptr;
888 }
889 ++t_row_base;
890 }
891 for (; rr < nb_base_functions; ++rr)
892 ++t_row_base;
893 ++t_coords;
894 ++t_w; // move to another integration weight
895 }
897};
898
899template <typename OpBase>
902 EntitiesFieldData::EntData &col_data) {
906 if (entsPtr) {
907 if (entsPtr->find(OpBase::getFEEntityHandle()) == entsPtr->end())
909 }
910 auto get_t_vec = [&](const int rr) {
911 std::array<double *, 3> ptrs;
912 for (auto i = 0; i != 3; ++i)
913 ptrs[i] = &OpBase::locMat(rr + i, i);
915 };
916 size_t nb_base_functions = row_data.getN().size2() / 3;
917 // // get element volume
918 const double vol = OpBase::getMeasure();
919 // get integration weights
920 auto t_w = OpBase::getFTensor0IntegrationWeight();
921 // get base function gradient on rows
922 auto t_row_base = row_data.getFTensor1N<3>();
923 // get coordinate at integration points
924 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
925 // loop over integration points
926 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
927 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
928 // take into account Jacobian
929 const double alpha = t_w * beta;
930 // loop over rows base functions
931 int rr = 0;
932 for (; rr != OpBase::nbRows / 3; rr++) {
933 // get column base functions gradient at gauss point gg
934 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
935 auto t_vec = get_t_vec(3 * rr);
936 // loop over columns
937 for (int cc = 0; cc != OpBase::nbCols / 3; cc++) {
938 // calculate element of local matrix
939 t_vec(i) += alpha * (t_row_base(k) * t_col_base(k));
940 ++t_col_base;
941 ++t_vec;
942 }
943 ++t_row_base;
944 }
945 for (; rr < nb_base_functions; ++rr)
946 ++t_row_base;
947 ++t_coords;
948 ++t_w; // move to another integration weight
949 }
951};
952
953template <int SPACE_DIM, int S, typename OpBase>
957 EntitiesFieldData::EntData &col_data) {
959
960 if (entsPtr) {
961 if (entsPtr->find(OpBase::getFEEntityHandle()) == entsPtr->end())
963 }
964
965 const size_t nb_row_dofs = row_data.getIndices().size();
966 const size_t nb_col_dofs = col_data.getIndices().size();
967
968 if (nb_row_dofs && nb_col_dofs) {
969
974
975 // get element volume
976 double vol = OpBase::getMeasure();
977
978 // get intergrayion weights
979 auto t_w = OpBase::getFTensor0IntegrationWeight();
980
981 // get derivatives of base functions on rows
982 auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
983
984 // Elastic stiffness tensor (4th rank tensor with minor and major
985 // symmetry)
986 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, S>(*matD);
987
988 // iterate over integration points
989 for (int gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
990
991 // calculate scalar weight times element volume
992 double a = t_w * vol;
993
994 // iterate over row base functions
995 int rr = 0;
996 for (; rr != OpBase::nbRows / SPACE_DIM; ++rr) {
997
998 // get sub matrix for the row
999 auto t_m = OpBase::template getLocMat<SPACE_DIM>(SPACE_DIM * rr);
1000
1002 // I mix up the indices here so that it behaves like a
1003 // Dg. That way I don't have to have a separate wrapper
1004 // class Christof_Expr, which simplifies things.
1005 t_rowD(l, j, k) = t_D(i, j, k, l) * (a * t_row_diff_base(i));
1006
1007 // get derivatives of base functions for columns
1008 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1009
1010 // iterate column base functions
1011 for (int cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1012
1013 // integrate block local stiffens matrix
1014 t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base(k);
1015
1016 // move to next column base function
1017 ++t_col_diff_base;
1018
1019 // move to next block of local stiffens matrix
1020 ++t_m;
1021 }
1022
1023 // move to next row base function
1024 ++t_row_diff_base;
1025 }
1026
1027 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
1028 ++t_row_diff_base;
1029
1030 // move to next integration weight
1031 ++t_w;
1032 ++t_D;
1033 }
1034 }
1035
1037}
1038
1039template <int SPACE_DIM, int S, typename OpBase>
1043 EntitiesFieldData::EntData &col_data) {
1045
1047
1052
1053 auto &row_hessian = row_data.getN(BaseDerivatives::SecondDerivative);
1054 auto &col_hessian = col_data.getN(BaseDerivatives::SecondDerivative);
1055
1056#ifndef NDEBUG
1057 if (row_hessian.size1() != OpBase::nbIntegrationPts) {
1058 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1059 "Wrong number of integration pts (%d != %d)",
1060 row_hessian.size1(), OpBase::nbIntegrationPts);
1061 }
1062 if (row_hessian.size2() !=
1064 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1065 "Wrong number of base functions (%d != %d)",
1066 row_hessian.size2() / (SPACE_DIM * SPACE_DIM),
1068 }
1069 if (row_hessian.size2() < OpBase::nbRows * SPACE_DIM * SPACE_DIM) {
1070 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1071 "Wrong number of base functions (%d < %d)", row_hessian.size2(),
1073 }
1074 if (col_hessian.size1() != OpBase::nbIntegrationPts) {
1075 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1076 "Wrong number of integration pts (%d != %d)",
1077 col_hessian.size1(), OpBase::nbIntegrationPts);
1078 }
1079 if (col_hessian.size2() < OpBase::nbCols * SPACE_DIM * SPACE_DIM) {
1080 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1081 "Wrong number of base functions (%d < %d)", col_hessian.size2(),
1083 }
1084#endif
1085
1086 // get element volume
1087 double vol = OpBase::getMeasure();
1088
1089 // get intergrayion weights
1090 auto t_w = OpBase::getFTensor0IntegrationWeight();
1091
1092 auto t_row_diff2 = getFTensor2SymmetricLowerFromPtr<SPACE_DIM>(
1093 &*row_hessian.data().begin());
1094
1095 // Elastic stiffness tensor (4th rank tensor with minor and major
1096 // symmetry)
1097 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, S>(*matD);
1098
1099 // iterate over integration points
1100 for (int gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1101
1102 // calculate scalar weight times element volume
1103 double a = t_w * vol;
1104
1105 // get sub matrix for the row
1106 auto m_ptr = &*OpBase::locMat.data().begin();
1107
1108 // iterate over row base functions
1109 int rr = 0;
1110 for (; rr != OpBase::nbRows; ++rr) {
1111
1113 t_rowD(k, l) = t_D(i, j, k, l) * (a * t_row_diff2(i, j));
1114
1115 // get derivatives of base functions for columns
1116 auto t_col_diff2 =
1117 getFTensor2SymmetricLowerFromPtr<SPACE_DIM>(&col_hessian(gg, 0));
1118
1119 // iterate column base functions
1120 for (int cc = 0; cc != OpBase::nbCols; ++cc) {
1121
1122 // integrate block local stiffens matrix
1123 *m_ptr += t_rowD(i, j) * t_col_diff2(i, j);
1124
1125 // move to next column base function
1126 ++t_col_diff2;
1127
1128 // move to next block of local stiffens matrix
1129 ++m_ptr;
1130 }
1131
1132 // move to next row base function
1133 ++t_row_diff2;
1134 }
1135
1136 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
1137 ++t_row_diff2;
1138
1139 // move to next integration weight
1140 ++t_w;
1141 ++t_D;
1142 }
1143 }
1144
1146}
1147
1148template <int SPACE_DIM, int S, typename OpBase>
1152 EntitiesFieldData::EntData &col_data) {
1154
1155 const size_t nb_row_dofs = row_data.getIndices().size();
1156 const size_t nb_col_dofs = col_data.getIndices().size();
1157
1158 if (nb_row_dofs && nb_col_dofs) {
1159
1164
1165 // get element volume
1166 double vol = OpBase::getMeasure();
1167
1168 // get intergrayion weights
1169 auto t_w = OpBase::getFTensor0IntegrationWeight();
1170
1171 // get derivatives of base functions on rows
1172 auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
1173
1174 // stiffness tensor (4th rank tensor)
1175 auto t_D =
1176 getFTensor4FromMat<SPACE_DIM, SPACE_DIM, SPACE_DIM, SPACE_DIM, S>(
1177 *matD);
1178
1179 // iterate over integration points
1180 for (int gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1181
1182 // calculate scalar weight times element volume
1183 double a = t_w * vol;
1184
1185 // iterate over row base functions
1186 int rr = 0;
1187 for (; rr != OpBase::nbRows / SPACE_DIM; ++rr) {
1188
1189 // get sub matrix for the row
1190 auto t_m = OpBase::template getLocMat<SPACE_DIM>(SPACE_DIM * rr);
1191
1192 // get derivatives of base functions for columns
1193 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1194
1195 // iterate column base functions
1196 for (int cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1197
1198 // integrate block local stiffens matrix
1199 t_m(i, k) +=
1200 a * (t_D(i, j, k, l) * (t_row_diff_base(j) * t_col_diff_base(l)));
1201
1202 // move to next column base function
1203 ++t_col_diff_base;
1204
1205 // move to next block of local stiffens matrix
1206 ++t_m;
1207 }
1208
1209 // move to next row base function
1210 ++t_row_diff_base;
1211 }
1212
1213 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
1214 ++t_row_diff_base;
1215
1216 // move to next integration weight
1217 ++t_w;
1218 ++t_D;
1219 }
1220 }
1221
1223}
1224
1225template <int SPACE_DIM, typename OpBase>
1228 EntitiesFieldData::EntData &col_data) {
1230
1231 auto t_w = this->getFTensor0IntegrationWeight();
1232
1233 size_t nb_base_functions = row_data.getN().size2() / 3;
1234 auto t_row_diff_base = row_data.getFTensor2DiffN<3, SPACE_DIM>();
1235
1236 const double alpha_constant = alphaConstant();
1237
1238 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1239 const double alpha = alpha_constant * this->getMeasure() * t_w;
1240
1241 size_t rr = 0;
1242 for (; rr != OpBase::nbRows; ++rr) {
1243 const double t_row_div_base = t_row_diff_base(i, i);
1244 auto t_col_base = col_data.getFTensor0N(gg, 0);
1245 for (size_t cc = 0; cc != OpBase::nbCols; ++cc) {
1246 this->locMat(rr, cc) += alpha * t_row_div_base * t_col_base;
1247 ++t_col_base;
1248 }
1249 ++t_row_diff_base;
1250 }
1251 for (; rr < nb_base_functions; ++rr)
1252 ++t_row_diff_base;
1253
1254 ++t_w;
1255 }
1256
1258}
1259
1260template <int SPACE_DIM, typename OpBase>
1263 EntitiesFieldData::EntData &col_data) {
1265
1266 auto t_w = this->getFTensor0IntegrationWeight();
1267
1268 size_t nb_base_functions = row_data.getN().size2() / 3;
1269 auto t_row_diff_base = row_data.getFTensor2DiffN<3, SPACE_DIM>();
1270 const double alpha_constant = alphaConstant();
1271 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1272
1273 const double alpha = alpha_constant * this->getMeasure() * t_w;
1274
1275 size_t rr = 0;
1276 for (; rr != OpBase::nbRows / SPACE_DIM; ++rr) {
1277 auto t_mat_diag = getFTensor1FromArrayDiag<SPACE_DIM, SPACE_DIM>(
1278 this->locMat, SPACE_DIM * rr);
1279 const double t_row_div_base = t_row_diff_base(i, i);
1280 auto t_col_base = col_data.getFTensor0N(gg, 0);
1281
1282 for (size_t cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1283 t_mat_diag(i) += alpha * t_row_div_base * t_col_base;
1284 ++t_col_base;
1285 ++t_mat_diag;
1286 }
1287
1288 ++t_row_diff_base;
1289 }
1290 for (; rr < nb_base_functions; ++rr)
1291 ++t_row_diff_base;
1292
1293 ++t_w;
1294 }
1295
1297}
1298
1299template <int SPACE_DIM, typename OpBase, CoordinateTypes COORDINATE_SYSTEM>
1303 EntitiesFieldData::EntData &col_data) {
1305
1306#ifndef NDEBUG
1307 if (OpBase::locMat.size2() % SPACE_DIM)
1308 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1309 "Number of rows in matrix should be multiple of space dimensions");
1310#endif
1311
1312 // When we move to C++17 add if constexpr()
1313 if constexpr (COORDINATE_SYSTEM == POLAR || COORDINATE_SYSTEM == SPHERICAL)
1314 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1315 "%s coordiante not implemented",
1316 CoordinateTypesNames[COORDINATE_SYSTEM]);
1317
1318 auto t_w = this->getFTensor0IntegrationWeight();
1319 auto t_coords = this->getFTensor1CoordsAtGaussPts();
1320 size_t nb_base_functions_row = row_data.getN().size2();
1321 auto t_row_base = row_data.getFTensor0N();
1322 const double vol = this->getMeasure();
1323 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1324
1325 const double alpha =
1326 alphaConstant(t_coords(0), t_coords(1), t_coords(2)) * t_w * vol;
1327
1328 size_t rr = 0;
1329 auto t_m = getFTensor1FromPtr<SPACE_DIM>(OpBase::locMat.data().data());
1330
1331 // When we move to C++17 add if constexpr()
1332 if constexpr (COORDINATE_SYSTEM == CARTESIAN) {
1333 for (; rr != OpBase::nbRows; ++rr) {
1334 const double r_val = alpha * t_row_base;
1335 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1336 for (size_t cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1337 t_m(i) += r_val * t_col_diff_base(i);
1338 ++t_col_diff_base;
1339 ++t_m;
1340 }
1341 ++t_row_base;
1342 }
1343 }
1344
1345 // When we move to C++17 add if constexpr()
1346 if constexpr (COORDINATE_SYSTEM == CYLINDRICAL) {
1347 for (; rr != OpBase::nbRows; ++rr) {
1348 const double r_val = alpha * t_row_base;
1349 auto t_col_base = col_data.getFTensor0N(gg, 0);
1350 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1351 for (size_t cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1352 t_m(i) += r_val * t_col_diff_base(i);
1353 t_m(0) += (r_val / t_coords(0)) * t_col_base;
1354 ++t_col_base;
1355 ++t_col_diff_base;
1356 ++t_m;
1357 }
1358 ++t_row_base;
1359 }
1360 }
1361
1362 for (; rr < nb_base_functions_row; ++rr)
1363 ++t_row_base;
1364
1365 ++t_w;
1366 ++t_coords;
1367 }
1368
1370}
1371
1372template <int SPACE_DIM, typename OpBase>
1376 EntitiesFieldData::EntData &col_data) {
1378
1379 auto t_w = this->getFTensor0IntegrationWeight();
1380
1381 size_t nb_base_functions = row_data.getN().size2() / 3;
1382 auto t_row_base = row_data.getFTensor1N<3>();
1383 auto &mat = this->locMat;
1384 const double alpha_constant = alphaConstant();
1385 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1386
1387 const double alpha = alpha_constant * this->getMeasure() * t_w;
1388
1389 size_t rr = 0;
1390 for (; rr != OpBase::nbRows; ++rr) {
1391 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1392 for (size_t cc = 0; cc != OpBase::nbCols; ++cc) {
1393 mat(rr, cc) += alpha * t_row_base(i) * t_col_diff_base(i);
1394 ++t_col_diff_base;
1395 }
1396 ++t_row_base;
1397 }
1398 for (; rr < nb_base_functions; ++rr)
1399 ++t_row_base;
1400
1401 ++t_w;
1402 }
1403
1405}
1406
1407template <int SPACE_DIM, typename OpBase>
1411 EntitiesFieldData::EntData &col_data) {
1413
1414 auto t_w = this->getFTensor0IntegrationWeight();
1415
1416 size_t nb_base_functions = row_data.getN().size2();
1417 auto t_row_base = row_data.getFTensor0N();
1418
1419 auto get_t_vec = [&](const int rr) {
1420 std::array<double *, SPACE_DIM> ptrs;
1421 for (auto i = 0; i != SPACE_DIM; ++i)
1422 ptrs[i] = &OpBase::locMat(rr + i, 0);
1424 };
1425
1426 const double alpha_constant = alphaConstant();
1427 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1428
1429 const double alpha = alpha_constant * this->getMeasure() * t_w;
1430
1431 size_t rr = 0;
1432 for (; rr != OpBase::nbRows / SPACE_DIM; ++rr) {
1433 auto t_vec = get_t_vec(SPACE_DIM * rr);
1434 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1435 for (size_t cc = 0; cc != OpBase::nbCols; ++cc) {
1436 t_vec(i) += alpha * t_row_base * t_col_diff_base(i);
1437 ++t_col_diff_base;
1438 ++t_vec;
1439 }
1440 ++t_row_base;
1441 }
1442 for (; rr < nb_base_functions; ++rr)
1443 ++t_row_base;
1444
1445 ++t_w;
1446 }
1447
1449}
1450
1451template <int SPACE_DIM, typename OpBase>
1454 EntitiesFieldData::EntData &col_data) {
1456
1457 auto t_w = this->getFTensor0IntegrationWeight();
1458
1459 size_t nb_base_functions = row_data.getN().size2() / 3;
1460 auto t_row_base = row_data.getFTensor1N<3>();
1461 const double alpha_constant = alphaConstant();
1462 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1463
1464 const double alpha = alpha_constant * this->getMeasure() * t_w;
1465
1466 size_t rr = 0;
1467 for (; rr != OpBase::nbRows / SPACE_DIM; ++rr) {
1468 auto t_mat_diag = getFTensor1FromArrayDiag<SPACE_DIM, SPACE_DIM>(
1469 this->locMat, SPACE_DIM * rr);
1470 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1471
1472 for (size_t cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1473 t_mat_diag(i) += alpha * t_row_base(j) * t_col_diff_base(j);
1474 ++t_col_diff_base;
1475 ++t_mat_diag;
1476 }
1477
1478 ++t_row_base;
1479 }
1480 for (; rr < nb_base_functions; ++rr)
1481 ++t_row_base;
1482
1483 ++t_w;
1484 }
1485
1487}
1488
1489template <int SPACE_DIM, typename OpBase>
1493 EntitiesFieldData::EntData &col_data) {
1495
1496 // get element volume
1497 const double vol = OpBase::getMeasure();
1498 // get integration weights
1499 auto t_w = OpBase::getFTensor0IntegrationWeight();
1500 // get base function gradient on rows
1501 auto t_row_base = row_data.getFTensor0N();
1502
1503 auto get_t_vec = [&](const int rr) {
1504 std::array<double *, SPACE_DIM> ptrs;
1505 for (auto i = 0; i != SPACE_DIM; ++i)
1506 ptrs[i] = &OpBase::locMat(rr, i);
1508 ptrs);
1509 };
1510
1511 auto t_grad_y = getFTensor1FromMat<SPACE_DIM>(*yGradPtr);
1513
1514 const double alpha_constant = alphaConstant();
1515 // loop over integration points
1516 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
1517 // take into account Jacobian
1518 const double alpha = t_w * vol * alpha_constant;
1519 // access local matrix
1520 auto t_vec = get_t_vec(0);
1521 // loop over rows base functions
1522 int rr = 0;
1523 for (; rr != OpBase::nbRows; rr++) {
1524 // get column base functions gradient at gauss point gg
1525 auto t_col_base = col_data.getFTensor0N(gg, 0);
1526 // loop over columns
1527 for (int cc = 0; cc != OpBase::nbCols / SPACE_DIM; cc++) {
1528 // calculate element of local matrix
1529 t_vec(i) += alpha * t_row_base * t_col_base * t_grad_y(i);
1530 ++t_col_base;
1531 ++t_vec;
1532 }
1533 ++t_row_base;
1534 }
1535 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
1536 ++t_row_base;
1537
1538 ++t_grad_y;
1539 ++t_w; // move to another integration weight
1540 }
1542};
1543
1544template <int SPACE_DIM, typename OpBase>
1548 EntitiesFieldData::EntData &col_data) {
1550
1551 // get element volume
1552 const double vol = OpBase::getMeasure();
1553 // get integration weights
1554 auto t_w = OpBase::getFTensor0IntegrationWeight();
1555 // get base function gradient on rows
1556 auto t_row_base = row_data.getFTensor0N();
1557
1558 auto t_u = getFTensor1FromMat<SPACE_DIM>(*uPtr);
1560 const double alpha_constant = alphaConstant();
1561 // loop over integration points
1562 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
1563 // take into account Jacobian
1564 const double alpha = t_w * vol * alpha_constant;
1565 // loop over rows base functions
1566 auto a_mat_ptr = &*OpBase::locMat.data().begin();
1567 int rr = 0;
1568 for (; rr != OpBase::nbRows; rr++) {
1569 // get column base functions gradient at gauss point gg
1570 auto t_diff_col_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1571 // loop over columns
1572 for (int cc = 0; cc != OpBase::nbCols; cc++) {
1573 // calculate element of local matrix
1574 (*a_mat_ptr) += alpha * t_row_base * t_diff_col_base(i) * t_u(i);
1575 ++t_diff_col_base;
1576 ++a_mat_ptr;
1577 }
1578 ++t_row_base;
1579 }
1580 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
1581 ++t_row_base;
1582
1583 ++t_u;
1584 ++t_w; // move to another integration weight
1585 }
1587};
1588
1589template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
1593 EntitiesFieldData::EntData &col_data) {
1595
1596 // get element volume
1597 const double vol = OpBase::getMeasure();
1598 // get integration weights
1599 auto t_w = OpBase::getFTensor0IntegrationWeight();
1600 // get base function gradient on rows
1601 auto t_row_base = row_data.getFTensor0N();
1602
1603 auto t_grad_y = getFTensor2FromMat<FIELD_DIM, SPACE_DIM>(*yGradPtr);
1606
1607 auto get_t_mat = [&](const int rr) {
1608 std::array<double *, FIELD_DIM * SPACE_DIM> ptrs;
1609 int s = 0;
1610 for (int j = 0; j != FIELD_DIM; ++j)
1611 for (auto i = 0; i != SPACE_DIM; ++i, ++s)
1612 ptrs[s] = &OpBase::locMat(rr + j, i);
1614 SPACE_DIM>(ptrs);
1615 };
1616
1617 const double alpha_constant = alphaConstant();
1618 // loop over integration points
1619 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
1620 // take into account Jacobian
1621 const double alpha = t_w * vol * alpha_constant;
1622
1623 // loop over rows base functions
1624 int rr = 0;
1625 for (; rr != OpBase::nbRows / FIELD_DIM; rr++) {
1626 // get column base functions gradient at gauss point gg
1627 auto t_col_base = col_data.getFTensor0N(gg, 0);
1628 // get mat
1629 auto t_mat = get_t_mat(FIELD_DIM * rr);
1630 // loop over columns
1631 for (int cc = 0; cc != OpBase::nbCols / SPACE_DIM; cc++) {
1632 // calculate element of local matrix
1633 t_mat(I, k) += (alpha * t_row_base * t_col_base) * t_grad_y(I, k);
1634 ++t_col_base;
1635 ++t_mat;
1636 }
1637 ++t_row_base;
1638 }
1639 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
1640 ++t_row_base;
1641
1642 ++t_grad_y;
1643 ++t_w; // move to another integration weight
1644 }
1646};
1647
1648template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
1652 EntitiesFieldData::EntData &col_data) {
1654
1655 // get element volume
1656 const double vol = OpBase::getMeasure();
1657 // get integration weights
1658 auto t_w = OpBase::getFTensor0IntegrationWeight();
1659 // get base function gradient on rows
1660 auto t_row_base = row_data.getFTensor0N();
1661
1662 auto get_t_mat = [&](const int rr) {
1663 std::array<double *, FIELD_DIM * FIELD_DIM> ptrs;
1664 int s = 0;
1665 for (int i = 0; i != FIELD_DIM; ++i)
1666 for (int j = 0; j != FIELD_DIM; ++j, ++s)
1667 ptrs[s] = &(OpBase::locMat(rr + i, j));
1669 FIELD_DIM>(ptrs);
1670 };
1671
1672 auto t_u = getFTensor1FromMat<SPACE_DIM>(*uPtr);
1674
1678
1679 const double alpha_constant = alphaConstant();
1680 // loop over integration points
1681 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
1682 // take into account Jacobian
1683 const double alpha = t_w * vol * alpha_constant;
1684
1685 // loop over rows base functions
1686 int rr = 0;
1687 for (; rr != OpBase::nbRows / FIELD_DIM; ++rr) {
1688 // get matrix vec
1689 auto t_mat = get_t_mat(FIELD_DIM * rr);
1690 // get column base functions gradient at gauss point gg
1691 auto t_diff_col_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1692 // loop over columns
1693 for (int cc = 0; cc != OpBase::nbCols / FIELD_DIM; ++cc) {
1694 t_mat(I, L) +=
1695 alpha * t_row_base * t_kd(I, L) * (t_diff_col_base(k) * t_u(k));
1696 ++t_mat;
1697 ++t_diff_col_base;
1698 }
1699 ++t_row_base;
1700 }
1701 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
1702 ++t_row_base;
1703
1704 ++t_u;
1705 ++t_w; // move to another integration weight
1706 }
1708};
1709
1710} // namespace MoFEM
1711
1712#endif //__BILINEAR_FORMS_INTEGRATORS_HPP__
static Index< 'I', 3 > I
constexpr double a
constexpr int SPACE_DIM
constexpr int FIELD_DIM
Kronecker Delta class symmetric.
static const char *const CoordinateTypesNames[]
Coordinate system names.
Definition: definitions.h:108
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
CoordinateTypes
Coodinate system.
Definition: definitions.h:114
@ CYLINDRICAL
Definition: definitions.h:117
@ POLAR
Definition: definitions.h:116
@ CARTESIAN
Definition: definitions.h:115
@ SPHERICAL
Definition: definitions.h:118
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
constexpr auto t_kd
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixScalarTimesDiv< SPACE_DIM, coord_type > OpMixScalarTimesDiv
IntegrationType
Form integrator integration types.
boost::function< double(const double, const double, const double)> ScalarFun
Scalar function type.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
constexpr int BASE_DIM
FTensor::Index< 'i', SPACE_DIM > i
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, 3 > OpGradGrad
Definition: heat_method.cpp:38
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
boost::function< double()> ConstantFun
Constant function type.
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
int rowSide
row side number
int colSide
column side number
int nbRows
number of dofs on rows
EntityType colType
column type
int nbIntegrationPts
number of integration points
MatrixDouble locMat
local entity block matrix
int nbCols
number if dof on column
int nbRowBaseFunctions
number or row base functions
EntityType rowType
row type
OpConvectiveTermLhsDuImpl(const std::string field_name_row, const std::string field_name_col, boost::shared_ptr< MatrixDouble > y_grad_ptr, ConstantFun alpha_fun=[]() { return 1;})
OpConvectiveTermLhsDuImpl(const std::string field_name_row, const std::string field_name_col, boost::shared_ptr< MatrixDouble > y_grad_ptr, ConstantFun alpha_fun=[]() { return 1;})
OpConvectiveTermLhsDyImpl(const std::string field_name_row, const std::string field_name_col, boost::shared_ptr< MatrixDouble > u_ptr, ConstantFun alpha_fun=[]() { return 1;})
OpConvectiveTermLhsDyImpl(const std::string field_name_row, const std::string field_name_col, boost::shared_ptr< MatrixDouble > u_ptr, ConstantFun alpha_fun=[]() { return 1;})
OpGradGradImpl(const std::string row_field_name, const std::string col_field_name, ScalarFun beta, boost::shared_ptr< Range > ents_ptr=nullptr)
OpGradGradImpl(const std::string row_field_name, const std::string col_field_name, ScalarFun beta, boost::shared_ptr< Range > ents_ptr=nullptr)
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)
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)
OpGradTensorGradImpl(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< MatrixDouble > mat_D)
OpMassImpl(const std::string row_field_name, const std::string col_field_name, ScalarFun beta, boost::shared_ptr< Range > ents_ptr=nullptr)
OpMassImpl(const std::string row_field_name, const std::string col_field_name, ScalarFun beta, boost::shared_ptr< Range > ents_ptr=nullptr)
OpMassImpl(const std::string row_field_name, const std::string col_field_name, ScalarFun beta, boost::shared_ptr< Range > ents_ptr=nullptr)
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)
OpMixDivTimesVecImpl(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)
OpMixScalarTimesDivImpl(const std::string row_field_name, const std::string col_field_name, ScalarFun alpha_fun, const bool assemble_transpose=false, const bool only_transpose=false)
OpMixTensorTimesGradImpl(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)
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)
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)