v0.14.0
Loading...
Searching...
No Matches
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#ifndef __BILINEAR_FORMS_INTEGRATORS_HPP__
14#define __BILINEAR_FORMS_INTEGRATORS_HPP__
15
16namespace MoFEM {
17
18template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, IntegrationType I,
19 typename OpBase>
21
22template <int SPACE_DIM, typename OpBase>
23struct OpGradGradImpl<1, 1, SPACE_DIM, GAUSS, OpBase> : public OpBase {
25 OpGradGradImpl(const std::string row_field_name,
26 const std::string col_field_name,
28 boost::shared_ptr<Range> ents_ptr = nullptr)
29 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
30 betaCoeff(beta) {
31 if (row_field_name == col_field_name)
32 this->sYmm = true;
33 }
34
35protected:
39};
40
41template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
44 OpGradGradImpl(const std::string row_field_name,
45 const std::string col_field_name,
47 boost::shared_ptr<Range> ents_ptr = nullptr)
48 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
49 betaCoeff(beta) {
50 if (row_field_name == col_field_name)
51 this->sYmm = true;
52 }
53
54protected:
58};
59
60template <int BASE_DIM, int FIELD_DIM, IntegrationType I, typename OpBase>
61struct OpMassImpl {};
62
63template <typename OpBase>
64struct OpMassImpl<1, 1, GAUSS, OpBase> : public OpBase {
65
67 const std::string row_field_name, const std::string col_field_name,
68 ScalarFun beta = [](double, double, double) constexpr { return 1; },
69 boost::shared_ptr<Range> ents_ptr = nullptr)
70 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
71 betaCoeff(beta) {
72 if (row_field_name == col_field_name)
73 this->sYmm = true;
74 }
75
76protected:
80};
81
82template <int FIELD_DIM, typename OpBase>
84 : public OpMassImpl<1, 1, GAUSS, OpBase> {
85 using OpMassImpl<1, 1, GAUSS, OpBase>::OpMassImpl;
86
87protected:
88 using OpMassImpl<1, 1, GAUSS, OpBase>::betaCoeff;
91};
92
93template <int FIELD_DIM, typename OpBase>
94struct OpMassImpl<3, FIELD_DIM, GAUSS, OpBase> : public OpBase {
95
96 OpMassImpl(const std::string row_field_name, const std::string col_field_name,
98 boost::shared_ptr<Range> ents_ptr = nullptr)
99 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
100 betaCoeff(beta) {
101 if (row_field_name == col_field_name)
102 this->sYmm = true;
103 }
104
105protected:
107 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
109};
110
111template <typename OpBase>
112struct OpMassImpl<3, 4, GAUSS, OpBase> : public OpBase {
113 OpMassImpl(const std::string row_field_name, const std::string col_field_name,
115 boost::shared_ptr<Range> ents_ptr = nullptr)
116 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
117 betaCoeff(beta) {
118 if (row_field_name == col_field_name)
119 this->sYmm = true;
120 }
121
122protected:
124 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
126};
127
128template <typename OpBase>
129struct OpMassImpl<3, 9, GAUSS, OpBase> : public OpBase {
130 OpMassImpl(const std::string row_field_name, const std::string col_field_name,
132 boost::shared_ptr<Range> ents_ptr = nullptr)
133 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
134 betaCoeff(beta) {
135 if (row_field_name == col_field_name)
136 this->sYmm = true;
137 }
138
139protected:
141 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
143};
144
145template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S, IntegrationType I,
146 typename OpBase>
148
149template <int SPACE_DIM, int S, typename OpBase>
151 : public OpBase {
153 OpGradSymTensorGradImpl(const std::string row_field_name,
154 const std::string col_field_name,
155 boost::shared_ptr<MatrixDouble> mat_D,
156 boost::shared_ptr<Range> ents_ptr = nullptr,
158 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
159 matD(mat_D), betaCoeff(beta) {
160 if (row_field_name == col_field_name)
161 this->sYmm = true;
162 }
163
164protected:
165 boost::shared_ptr<MatrixDouble> matD;
167 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
169};
170
171template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S, IntegrationType I,
172 typename OpBase>
174
175template <int SPACE_DIM, int S, typename OpBase>
177 : public OpBase {
179 OpGradTensorGradImpl(const std::string row_field_name,
180 const std::string col_field_name,
181 boost::shared_ptr<MatrixDouble> mat_D,
183 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL), matD(mat_D),
184 betaCoeff(beta) {}
185
186protected:
187 boost::shared_ptr<MatrixDouble> matD;
189 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
191};
192
193template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S, IntegrationType I,
194 typename OpBase>
196
197template <int SPACE_DIM, int S, typename OpBase>
199 : public OpBase {
201 OpGradGradSymTensorGradGradImpl(const std::string row_field_name,
202 const std::string col_field_name,
203 boost::shared_ptr<MatrixDouble> mat_D,
204 boost::shared_ptr<Range> ents_ptr = nullptr)
205 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
206 matD(mat_D) {
207 if (row_field_name == col_field_name)
208 this->sYmm = true;
209 }
210
211protected:
212 boost::shared_ptr<MatrixDouble> matD;
213 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
215};
216
217template <int SPACE_DIM, IntegrationType I, typename OpBase>
219
220template <int SPACE_DIM, typename OpBase>
222 OpMixDivTimesScalarImpl(const std::string row_field_name,
223 const std::string col_field_name,
224 ConstantFun alpha_fun,
225 const bool assemble_transpose = false,
226 const bool only_transpose = false)
227 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
228 alphaConstant(alpha_fun) {
229 this->assembleTranspose = assemble_transpose;
230 this->onlyTranspose = only_transpose;
231 }
232
233protected:
236 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
238};
239
240template <int SPACE_DIM, IntegrationType I, typename OpBase,
241 CoordinateTypes CoordSys>
243
244template <int SPACE_DIM, typename OpBase, CoordinateTypes CoordSys>
246 : public OpBase {
247
248 OpMixDivTimesVecImpl(const std::string row_field_name,
249 const std::string col_field_name, ConstantFun alpha_fun,
250 const bool assemble_transpose,
251 const bool only_transpose = false)
252 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
253 alphaConstant(alpha_fun) {
254 this->assembleTranspose = assemble_transpose;
255 this->onlyTranspose = only_transpose;
256 }
257
258 OpMixDivTimesVecImpl(const std::string row_field_name,
259 const std::string col_field_name, ConstantFun alpha_fun,
260 ScalarFun beta_fun, const bool assemble_transpose,
261 const bool only_transpose = false)
262 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
263 alphaConstant(alpha_fun), betaCoeff(beta_fun) {
264 this->assembleTranspose = assemble_transpose;
265 this->onlyTranspose = only_transpose;
266 }
267
268protected:
270
271 ConstantFun alphaConstant = []() constexpr { return 1; };
272 ScalarFun betaCoeff = [](double, double, double) constexpr { return 1; };
273
274 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
276};
277
278template <int SPACE_DIM, IntegrationType I, typename OpBase,
279 CoordinateTypes COORDINATE_SYSTEM>
281
282template <int SPACE_DIM, typename OpBase, CoordinateTypes COORDINATE_SYSTEM>
283struct OpMixScalarTimesDivImpl<SPACE_DIM, GAUSS, OpBase, COORDINATE_SYSTEM>
284 : public OpBase {
286 const std::string row_field_name, const std::string col_field_name,
287 ScalarFun alpha_fun = [](double, double, double) constexpr { return 1; },
288 const bool assemble_transpose = false, const bool only_transpose = false)
289 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
290 alphaConstant(alpha_fun) {
291 this->assembleTranspose = assemble_transpose;
292 this->onlyTranspose = only_transpose;
293 this->sYmm = false;
294 }
295
296protected:
299 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
301};
302
303template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, IntegrationType I,
304 typename OpBase>
306
307template <int SPACE_DIM, typename OpBase>
309 : public OpBase {
310 OpMixVectorTimesGradImpl(const std::string row_field_name,
311 const std::string col_field_name,
312 ConstantFun alpha_fun,
313 const bool assemble_transpose = false,
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
321protected:
324 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
326};
327
328template <int SPACE_DIM, typename OpBase>
330 : public OpBase {
331 OpMixVectorTimesGradImpl(const std::string row_field_name,
332 const std::string col_field_name,
333 ConstantFun alpha_fun,
334 const bool assemble_transpose = false,
335 const bool only_transpose = false)
336 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
337 alphaConstant(alpha_fun) {
338 this->assembleTranspose = assemble_transpose;
339 this->onlyTranspose = only_transpose;
340 }
341
342protected:
345 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
347};
348
349template <int SPACE_DIM, IntegrationType I, typename OpBase>
351
352template <int SPACE_DIM, typename OpBase>
354
355 OpMixTensorTimesGradImpl(const std::string row_field_name,
356 const std::string col_field_name,
357 ConstantFun alpha_fun, const bool assemble_transpose,
358 const bool only_transpose = false)
359 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
360 alphaConstant(alpha_fun) {
361 this->assembleTranspose = assemble_transpose;
362 this->onlyTranspose = only_transpose;
363 }
364
365 OpMixTensorTimesGradImpl(const std::string row_field_name,
366 const std::string col_field_name,
367 ConstantFun alpha_fun, ScalarFun beta_coeff,
368 const bool assemble_transpose,
369 const bool only_transpose = false)
370 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
371 alphaConstant(alpha_fun), betaCoeff(beta_coeff) {
372 this->assembleTranspose = assemble_transpose;
373 this->onlyTranspose = only_transpose;
374 }
375
376protected:
379 ConstantFun alphaConstant = []() constexpr { return 1; };
380 ScalarFun betaCoeff = [](double, double, double) constexpr { return 1; };
381 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
383};
384
385/**
386 * @brief Bilinear integrator form
387 * @ingroup mofem_forms
388 *
389 * @tparam EleOp
390 * @tparam A
391 * @tparam I
392 */
393template <typename EleOp>
394template <AssemblyType A>
395template <IntegrationType I>
397
398 /**
399 * @brief Integrate \f$(v_{,i},\beta(\mathbf{x}) u_{,j}))_\Omega\f$
400 * @ingroup mofem_forms
401 *
402 * @tparam SPACE_DIM
403 */
404 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM>
406
407 /**
408 * @brief Integrate \f$(v_i,\beta(\mathbf{x}) u_j)_\Omega\f$
409 * @ingroup mofem_forms
410 *
411 * @note That FIELD_DIM = 4 or 9 is assumed that OpMass is for tensorial field
412 * 2x2 or 3x3
413 *
414 * @todo It should be considered another template parameter RANK which will
415 * allow to distinguish between scalar, vectorial and tensorial fields
416 *
417 * @tparam BASE_DIM dimension of base
418 * @tparam FIELD_DIM dimension of field
419 */
420 template <int BASE_DIM, int FIELD_DIM>
422
423 /**
424 * @brief Integrate \f$(v_k,D_{ijkl} u_{,l})_\Omega\f$
425 *
426 * \note \f$D_{ijkl}\f$ is * tensor with minor & major symmetry
427 *
428 * @ingroup mofem_forms
429 *
430 * @tparam SPACE_DIM
431 * @tparam S
432 */
433 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S = 0>
436
437 /**
438 * @brief Integrate \f$(v_{,ij},D_{ijkl} u_{,kl})_\Omega\f$
439 *
440 * \note \f$D_{ijkl}\f$ is * tensor with no symmetries
441 *
442 * @ingroup mofem_forms
443 *
444 * @tparam SPACE_DIM
445 * @tparam S
446 */
447 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S = 0>
450 OpBase>;
451
452 /**
453 * @brief Integrate \f$(v_{,j},D_{ijkl} u_{,l})_\Omega\f$
454 *
455 * \note \f$D_{ijkl}\f$ is * tensor with no symmetries
456 *
457 * @ingroup mofem_forms
458 *
459 * @tparam SPACE_DIM
460 * @tparam S
461 */
462 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S = 0>
465
466 /**
467 * @brief Integrate \f$(\lambda_{i,i},u)_\Omega\f$
468 *
469 * @tparam SPACE_DIM
470 */
471 template <int SPACE_DIM>
473
474 /**
475 * @brief Integrate \f$(\lambda_{ij,j},u_{i})_\Omega\f$
476 *
477 * @tparam SPACE_DIM
478 */
479 template <int SPACE_DIM, CoordinateTypes CoordSys = CARTESIAN>
481
482 /**
483 * @brief Integrate \f$(\lambda,u_{i,i})_\Omega\f$
484 *
485 * @tparam SPACE_DIM
486 */
487 template <int SPACE_DIM, CoordinateTypes COORDINATE_SYSTEM = CARTESIAN>
490
491 /**
492 * @brief Integrate \f$(\lambda_{i},u_{,j})_\Omega\f$
493 *
494 * @tparam SPACE_DIM
495 */
496 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM>
499
500 /**
501 * @brief Integrate \f$(\lambda_{ij},u_{i,j})_\Omega\f$
502 *
503 * @tparam SPACE_DIM
504 */
505 template <int SPACE_DIM>
507};
508
509template <int SPACE_DIM, typename OpBase>
512 EntitiesFieldData::EntData &col_data) {
514
515 // get element volume
516 const double vol = OpBase::getMeasure();
517 // get integration weights
518 auto t_w = OpBase::getFTensor0IntegrationWeight();
519 // get base function gradient on rows
520 auto t_row_grad = row_data.getFTensor1DiffN<SPACE_DIM>();
521 // get coordinate at integration points
522 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
523
524 // loop over integration points
525 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
526 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
527 // take into account Jacobian
528 const double alpha = t_w * beta;
529 // loop over ros base functions
530 int rr = 0;
531 for (; rr != OpBase::nbRows; rr++) {
532 // get column base functions gradient at gauss point gg
533 auto t_col_grad = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
534 // loop over columns
535 for (int cc = 0; cc != OpBase::nbCols; cc++) {
536 // calculate element of local matrix
537 OpBase::locMat(rr, cc) += alpha * (t_row_grad(i) * t_col_grad(i));
538 ++t_col_grad; // move to another gradient of base function
539 // on column
540 }
541 ++t_row_grad; // move to another element of gradient of base
542 // function on row
543 }
544 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
545 ++t_row_grad;
546
547 ++t_coords;
548 ++t_w; // move to another integration weight
549 }
551}
552
553template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
557 EntitiesFieldData::EntData &col_data) {
559
560 // get element volume
561 const double vol = OpBase::getMeasure();
562 // get integration weights
563 auto t_w = OpBase::getFTensor0IntegrationWeight();
564 // get base function gradient on rows
565 auto t_row_grad = row_data.getFTensor1DiffN<SPACE_DIM>();
566 // get coordinate at integration points
567 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
568
571
572 auto get_t_vec = [&](const int rr) {
573 std::array<double *, FIELD_DIM> ptrs;
574 for (auto i = 0; i != FIELD_DIM; ++i)
575 ptrs[i] = &OpBase::locMat(rr + i, i);
577 ptrs);
578 };
579
580 // loop over integration points
581 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
582 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
583 // take into account Jacobian
584 const double alpha = t_w * beta;
585 // loop over ros base functions
586 int rr = 0;
587 for (; rr != OpBase::nbRows / FIELD_DIM; rr++) {
588 // get diag vec
589 auto t_vec = get_t_vec(rr * FIELD_DIM);
590 // get column base functions gradient at gauss point gg
591 auto t_col_grad = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
592 // loop over columns
593 for (int cc = 0; cc != OpBase::nbCols / FIELD_DIM; cc++) {
594 // calculate element of local matrix
595 t_vec(j) += alpha * (t_row_grad(i) * t_col_grad(i));
596 ++t_col_grad; // move to another gradient of base function
597 // on column
598 ++t_vec;
599 }
600 ++t_row_grad; // move to another element of gradient of base
601 // function on row
602 }
603 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
604 ++t_row_grad;
605
606 ++t_coords;
607 ++t_w; // move to another integration weight
608 }
610}
611
612template <typename OpBase>
615 EntitiesFieldData::EntData &col_data) {
617
618#ifndef NDEBUG
619 auto log_error = [&]() {
620 MOFEM_LOG("SELF", Sev::error) << "Row side " << OpBase::rowSide << " "
621 << CN::EntityTypeName(OpBase::rowType);
622 MOFEM_LOG("SELF", Sev::error) << "Col side " << OpBase::colSide << " "
623 << CN::EntityTypeName(OpBase::colType);
624 };
625
626 if (row_data.getN().size2() < OpBase::nbRows) {
627 log_error();
628 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
629 "Wrong number of base functions on rows %d < %d",
630 row_data.getN().size2(), OpBase::nbRows);
631 }
632 if (col_data.getN().size2() < OpBase::nbCols) {
633 log_error();
634 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
635 "Wrong number of base functions on cols %d < %d",
636 col_data.getN().size2(), OpBase::nbCols);
637 }
638 if (row_data.getN().size1() != OpBase::nbIntegrationPts) {
639 log_error();
640 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
641 "Wrong number of integration points on rows %d != %d",
642 row_data.getN().size1(), OpBase::nbIntegrationPts);
643 }
644 if (col_data.getN().size1() != OpBase::nbIntegrationPts) {
645 log_error();
646 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
647 "Wrong number of integration points on cols %d != %d",
648 col_data.getN().size1(), OpBase::nbIntegrationPts);
649 }
650#endif
651
652 // get element volume
653 const double vol = OpBase::getMeasure();
654 // get integration weights
655 auto t_w = OpBase::getFTensor0IntegrationWeight();
656 // get base function gradient on rows
657 auto t_row_base = row_data.getFTensor0N();
658 // get coordinate at integration points
659 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
660 // loop over integration points
661 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
662 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
663 // take into account Jacobian
664 const double alpha = t_w * beta;
665 // loop over rows base functions
666 auto a_mat_ptr = &*OpBase::locMat.data().begin();
667 int rr = 0;
668 for (; rr != OpBase::nbRows; rr++) {
669 // get column base functions gradient at gauss point gg
670 auto t_col_base = col_data.getFTensor0N(gg, 0);
671 // loop over columns
672 for (int cc = 0; cc != OpBase::nbCols; cc++) {
673 // calculate element of local matrix
674 *a_mat_ptr += alpha * (t_row_base * t_col_base);
675 ++t_col_base;
676 ++a_mat_ptr;
677 }
678 ++t_row_base;
679 }
680 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
681 ++t_row_base;
682 ++t_coords;
683 ++t_w; // move to another integration weight
684 }
686};
687
688template <int FIELD_DIM, typename OpBase>
691 EntitiesFieldData::EntData &col_data) {
693 // get element volume
694 const double vol = OpBase::getMeasure();
695 // get integration weights
696 auto t_w = OpBase::getFTensor0IntegrationWeight();
697 // get base function gradient on rows
698 auto t_row_base = row_data.getFTensor0N();
699 // get coordinate at integration points
700 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
701
703 auto get_t_vec = [&](const int rr) {
704 std::array<double *, FIELD_DIM> ptrs;
705 for (auto i = 0; i != FIELD_DIM; ++i)
706 ptrs[i] = &OpBase::locMat(rr + i, i);
708 ptrs);
709 };
710
711 // loop over integration points
712 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
713 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
714 // take into account Jacobian
715 const double alpha = t_w * beta;
716 // loop over rows base functions
717 int rr = 0;
718 for (; rr != OpBase::nbRows / FIELD_DIM; rr++) {
719 // get column base functions gradient at gauss point gg
720 auto t_col_base = col_data.getFTensor0N(gg, 0);
721 // get mat vec
722 auto t_vec = get_t_vec(FIELD_DIM * rr);
723 // loop over columns
724 for (int cc = 0; cc != OpBase::nbCols / FIELD_DIM; cc++) {
725 // calculate element of local matrix
726 t_vec(i) += alpha * (t_row_base * t_col_base);
727 ++t_col_base;
728 ++t_vec;
729 }
730 ++t_row_base;
731 }
732 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
733 ++t_row_base;
734 ++t_coords;
735 ++t_w; // move to another integration weight
736 }
738};
739
740template <int FIELD_DIM, typename OpBase>
743 EntitiesFieldData::EntData &col_data) {
746 size_t nb_base_functions = row_data.getN().size2() / 3;
747 // // get element volume
748 const double vol = OpBase::getMeasure();
749 // get integration weights
750 auto t_w = OpBase::getFTensor0IntegrationWeight();
751 // get base function gradient on rows
752 auto t_row_base = row_data.getFTensor1N<3>();
753 // get coordinate at integration points
754 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
755 // loop over integration points
756 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
757 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
758 // take into account Jacobian
759 const double alpha = t_w * beta;
760 // loop over rows base functions
761 auto a_mat_ptr = &*OpBase::locMat.data().begin();
762 int rr = 0;
763 for (; rr != OpBase::nbRows; rr++) {
764 // get column base functions gradient at gauss point gg
765 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
766 // loop over columns
767 for (int cc = 0; cc != OpBase::nbCols; cc++) {
768 // calculate element of local matrix
769 (*a_mat_ptr) += alpha * (t_row_base(i) * t_col_base(i));
770 ++t_col_base;
771 ++a_mat_ptr;
772 }
773 ++t_row_base;
774 }
775 for (; rr < nb_base_functions; ++rr)
776 ++t_row_base;
777 ++t_coords;
778 ++t_w; // move to another integration weight
779 }
781};
782
783template <typename OpBase>
786 EntitiesFieldData::EntData &col_data) {
790 auto get_t_vec = [&](const int rr) {
792 &OpBase::locMat(rr + 0, 0), &OpBase::locMat(rr + 1, 1)};
793 };
794 size_t nb_base_functions = row_data.getN().size2() / 3;
795 // // get element volume
796 const double vol = OpBase::getMeasure();
797 // get integration weights
798 auto t_w = OpBase::getFTensor0IntegrationWeight();
799 // get base function gradient on rows
800 auto t_row_base = row_data.getFTensor1N<3>();
801 // get coordinate at integration points
802 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
803 // loop over integration points
804 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
805 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
806 // take into account Jacobian
807 const double alpha = t_w * beta;
808 // loop over rows base functions
809 int rr = 0;
810 for (; rr != OpBase::nbRows / 2; rr++) {
811 // get column base functions gradient at gauss point gg
812 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
813 auto t_vec = get_t_vec(2 * rr);
814 // loop over columns
815 for (int cc = 0; cc != OpBase::nbCols / 2; cc++) {
816 // calculate element of local matrix
817 t_vec(I) += alpha * (t_row_base(k) * t_col_base(k));
818 ++t_col_base;
819 ++t_vec;
820 }
821 ++t_row_base;
822 }
823 for (; rr < nb_base_functions; ++rr)
824 ++t_row_base;
825 ++t_coords;
826 ++t_w; // move to another integration weight
827 }
829}
830
831template <typename OpBase>
834 EntitiesFieldData::EntData &col_data) {
838 auto get_t_vec = [&](const int rr) {
840 &OpBase::locMat(rr + 0, 0), &OpBase::locMat(rr + 1, 1),
841 &OpBase::locMat(rr + 2, 2)};
842 };
843 size_t nb_base_functions = row_data.getN().size2() / 3;
844 // // get element volume
845 const double vol = OpBase::getMeasure();
846 // get integration weights
847 auto t_w = OpBase::getFTensor0IntegrationWeight();
848 // get base function gradient on rows
849 auto t_row_base = row_data.getFTensor1N<3>();
850 // get coordinate at integration points
851 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
852 // loop over integration points
853 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
854 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
855 // take into account Jacobian
856 const double alpha = t_w * beta;
857 // loop over rows base functions
858 int rr = 0;
859 for (; rr != OpBase::nbRows / 3; rr++) {
860 // get column base functions gradient at gauss point gg
861 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
862 auto t_vec = get_t_vec(3 * rr);
863 // loop over columns
864 for (int cc = 0; cc != OpBase::nbCols / 3; cc++) {
865 // calculate element of local matrix
866 t_vec(i) += alpha * (t_row_base(k) * t_col_base(k));
867 ++t_col_base;
868 ++t_vec;
869 }
870 ++t_row_base;
871 }
872 for (; rr < nb_base_functions; ++rr)
873 ++t_row_base;
874 ++t_coords;
875 ++t_w; // move to another integration weight
876 }
878};
879
880template <int SPACE_DIM, int S, typename OpBase>
884 EntitiesFieldData::EntData &col_data) {
886
887 const size_t nb_row_dofs = row_data.getIndices().size();
888 const size_t nb_col_dofs = col_data.getIndices().size();
889
890 if (nb_row_dofs && nb_col_dofs) {
891
892 const bool diag = (row_data.getFieldEntities()[0]->getLocalUniqueId() ==
893 col_data.getFieldEntities()[0]->getLocalUniqueId());
894
899
900 // get element volume
901 double vol = OpBase::getMeasure();
902
903 // get intergrayion weights
904 auto t_w = OpBase::getFTensor0IntegrationWeight();
905
906 // get derivatives of base functions on rows
907 auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
908
909 // Elastic stiffness tensor (4th rank tensor with minor and major
910 // symmetry)
911 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, S>(*matD);
912
913 // get coordinate at integration points
914 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
915
916 // iterate over integration points
917 for (int gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
918
919 // calculate scalar weight times element volume
920 double a = t_w * vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
921
922 // iterate over row base functions
923 int rr = 0;
924 for (; rr != OpBase::nbRows / SPACE_DIM; ++rr) {
925
926 // get sub matrix for the row
927 auto t_m = OpBase::template getLocMat<SPACE_DIM>(SPACE_DIM * rr);
928
930 // I mix up the indices here so that it behaves like a
931 // Dg. That way I don't have to have a separate wrapper
932 // class Christof_Expr, which simplifies things.
933 t_rowD(l, j, k) = t_D(i, j, k, l) * (a * t_row_diff_base(i));
934
935 // get derivatives of base functions for columns
936 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
937
938 // iterate column base functions
939 const auto nb_cols = (diag) ? rr : OpBase::nbCols / SPACE_DIM - 1;
940 for (int cc = 0; cc <= nb_cols; ++cc) {
941
942 // integrate block local stiffens matrix
943 t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base(k);
944
945 // move to next column base function
946 ++t_col_diff_base;
947
948 // move to next block of local stiffens matrix
949 ++t_m;
950 }
951
952 // move to next row base function
953 ++t_row_diff_base;
954 }
955
956 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
957 ++t_row_diff_base;
958
959 // move to next integration weight
960 ++t_w;
961 ++t_D;
962 ++t_coords;
963 }
964
965 // Copy symmetry
966 if (diag) {
967 for (int rr = 0; rr != OpBase::nbRows / SPACE_DIM - 1; ++rr) {
968 auto t_m_rr = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
969 this->locMat, SPACE_DIM * rr, SPACE_DIM * (rr + 1));
970 auto t_m_cc = getFTensor2FromArray<SPACE_DIM, SPACE_DIM>(
971 this->locMat, SPACE_DIM * (rr + 1), SPACE_DIM * rr,
973 for (int cc = rr + 1; cc < OpBase::nbCols / SPACE_DIM; ++cc) {
974 t_m_rr(i, j) = t_m_cc(j, i);
975 ++t_m_rr;
976 ++t_m_cc;
977 }
978 }
979 }
980 }
981
983}
984
985template <int SPACE_DIM, int S, typename OpBase>
989 EntitiesFieldData::EntData &col_data) {
991
993
998
999 auto &row_hessian = row_data.getN(BaseDerivatives::SecondDerivative);
1000 auto &col_hessian = col_data.getN(BaseDerivatives::SecondDerivative);
1001
1002#ifndef NDEBUG
1003 if (row_hessian.size1() != OpBase::nbIntegrationPts) {
1004 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1005 "Wrong number of integration pts (%d != %d)",
1006 row_hessian.size1(), OpBase::nbIntegrationPts);
1007 }
1008 if (row_hessian.size2() !=
1010 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1011 "Wrong number of base functions (%d != %d)",
1012 row_hessian.size2() / (SPACE_DIM * SPACE_DIM),
1014 }
1015 if (row_hessian.size2() < OpBase::nbRows * SPACE_DIM * SPACE_DIM) {
1016 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1017 "Wrong number of base functions (%d < %d)", row_hessian.size2(),
1019 }
1020 if (col_hessian.size1() != OpBase::nbIntegrationPts) {
1021 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1022 "Wrong number of integration pts (%d != %d)",
1023 col_hessian.size1(), OpBase::nbIntegrationPts);
1024 }
1025 if (col_hessian.size2() < OpBase::nbCols * SPACE_DIM * SPACE_DIM) {
1026 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1027 "Wrong number of base functions (%d < %d)", col_hessian.size2(),
1029 }
1030#endif
1031
1032 // get element volume
1033 double vol = OpBase::getMeasure();
1034
1035 // get intergrayion weights
1036 auto t_w = OpBase::getFTensor0IntegrationWeight();
1037
1038 auto t_row_diff2 = getFTensor2SymmetricLowerFromPtr<SPACE_DIM>(
1039 &*row_hessian.data().begin());
1040
1041 // Elastic stiffness tensor (4th rank tensor with minor and major
1042 // symmetry)
1043 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, S>(*matD);
1044
1045 // iterate over integration points
1046 for (int gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1047
1048 // calculate scalar weight times element volume
1049 double a = t_w * vol;
1050
1051 // get sub matrix for the row
1052 auto m_ptr = &*OpBase::locMat.data().begin();
1053
1054 // iterate over row base functions
1055 int rr = 0;
1056 for (; rr != OpBase::nbRows; ++rr) {
1057
1059 t_rowD(k, l) = t_D(i, j, k, l) * (a * t_row_diff2(i, j));
1060
1061 // get derivatives of base functions for columns
1062 auto t_col_diff2 =
1063 getFTensor2SymmetricLowerFromPtr<SPACE_DIM>(&col_hessian(gg, 0));
1064
1065 // iterate column base functions
1066 for (int cc = 0; cc != OpBase::nbCols; ++cc) {
1067
1068 // integrate block local stiffens matrix
1069 *m_ptr += t_rowD(i, j) * t_col_diff2(i, j);
1070
1071 // move to next column base function
1072 ++t_col_diff2;
1073
1074 // move to next block of local stiffens matrix
1075 ++m_ptr;
1076 }
1077
1078 // move to next row base function
1079 ++t_row_diff2;
1080 }
1081
1082 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
1083 ++t_row_diff2;
1084
1085 // move to next integration weight
1086 ++t_w;
1087 ++t_D;
1088 }
1089 }
1090
1092}
1093
1094template <int SPACE_DIM, int S, typename OpBase>
1098 EntitiesFieldData::EntData &col_data) {
1100
1101 const size_t nb_row_dofs = row_data.getIndices().size();
1102 const size_t nb_col_dofs = col_data.getIndices().size();
1103
1104 if (nb_row_dofs && nb_col_dofs) {
1105
1110
1111 // get element volume
1112 double vol = OpBase::getMeasure();
1113
1114 // get intergrayion weights
1115 auto t_w = OpBase::getFTensor0IntegrationWeight();
1116
1117 // get derivatives of base functions on rows
1118 auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
1119
1120 // stiffness tensor (4th rank tensor)
1121 auto t_D =
1122 getFTensor4FromMat<SPACE_DIM, SPACE_DIM, SPACE_DIM, SPACE_DIM, S>(
1123 *matD);
1124
1125 // get coordinate at integration points
1126 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
1127
1128 // iterate over integration points
1129 for (int gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1130
1131 // calculate scalar weight times element volume
1132 double a = t_w * vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
1133
1134 // iterate over row base functions
1135 int rr = 0;
1136 for (; rr != OpBase::nbRows / SPACE_DIM; ++rr) {
1137
1138 // get sub matrix for the row
1139 auto t_m = OpBase::template getLocMat<SPACE_DIM>(SPACE_DIM * rr);
1140
1141 // calculate row
1143 t_row(i, k, l) = t_D(i, j, k, l) * (a * t_row_diff_base(j));
1144
1145 // get derivatives of base functions for columns
1146 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1147
1148 // iterate column base functions
1149 for (int cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1150
1151 // integrate block local stiffens matrix
1152 t_m(i, k) += t_row(i, k, l) * t_col_diff_base(l);
1153
1154 // move to next column base function
1155 ++t_col_diff_base;
1156
1157 // move to next block of local stiffens matrix
1158 ++t_m;
1159 }
1160
1161 // move to next row base function
1162 ++t_row_diff_base;
1163 }
1164
1165 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
1166 ++t_row_diff_base;
1167
1168 // move to next integration weight
1169 ++t_w;
1170 ++t_D;
1171 ++t_coords;
1172 }
1173 }
1174
1176}
1177
1178template <int SPACE_DIM, typename OpBase>
1181 EntitiesFieldData::EntData &col_data) {
1183
1184 auto t_w = this->getFTensor0IntegrationWeight();
1185
1186 size_t nb_base_functions = row_data.getN().size2() / 3;
1187 auto t_row_diff_base = row_data.getFTensor2DiffN<3, SPACE_DIM>();
1188
1189 const double alpha_constant = alphaConstant();
1190
1191 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1192 const double alpha = alpha_constant * this->getMeasure() * t_w;
1193
1194 size_t rr = 0;
1195 for (; rr != OpBase::nbRows; ++rr) {
1196 const double t_row_div_base = t_row_diff_base(i, i);
1197 auto t_col_base = col_data.getFTensor0N(gg, 0);
1198 for (size_t cc = 0; cc != OpBase::nbCols; ++cc) {
1199 this->locMat(rr, cc) += alpha * t_row_div_base * t_col_base;
1200 ++t_col_base;
1201 }
1202 ++t_row_diff_base;
1203 }
1204 for (; rr < nb_base_functions; ++rr)
1205 ++t_row_diff_base;
1206
1207 ++t_w;
1208 }
1209
1211}
1212
1213template <int SPACE_DIM, typename OpBase, CoordinateTypes CoordSys>
1217 EntitiesFieldData::EntData &col_data) {
1219
1220 auto t_w = this->getFTensor0IntegrationWeight();
1221 auto t_coords = this->getFTensor1CoordsAtGaussPts();
1222
1223 size_t nb_base_functions = row_data.getN().size2() / 3;
1224 auto t_row_diff_base = row_data.getFTensor2DiffN<3, SPACE_DIM>();
1225 auto t_row_base = row_data.getFTensor1N<3>();
1226 const double alpha_constant = alphaConstant() * this->getMeasure();
1227 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1228
1229 const double alpha =
1230 alpha_constant * t_w * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
1231
1232 size_t rr = 0;
1233 for (; rr != OpBase::nbRows / SPACE_DIM; ++rr) {
1234 auto t_mat_diag = getFTensor1FromArrayDiag<SPACE_DIM, SPACE_DIM>(
1235 this->locMat, SPACE_DIM * rr);
1236 const double t_row_div_base = t_row_diff_base(i, i);
1237 auto t_col_base = col_data.getFTensor0N(gg, 0);
1238
1239 for (size_t cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1240 t_mat_diag(i) += alpha * t_row_div_base * t_col_base;
1241 if constexpr (CoordSys == CYLINDRICAL) {
1242 t_mat_diag(i) += alpha * (t_row_base(0) / t_coords(0)) * t_col_base;
1243 }
1244 ++t_col_base;
1245 ++t_mat_diag;
1246 }
1247 ++t_row_base;
1248 ++t_row_diff_base;
1249 }
1250 for (; rr < nb_base_functions; ++rr) {
1251 ++t_row_diff_base;
1252 ++t_row_base;
1253 }
1254
1255 ++t_w;
1256 ++t_coords;
1257 }
1258
1260}
1261
1262template <int SPACE_DIM, typename OpBase, CoordinateTypes COORDINATE_SYSTEM>
1266 EntitiesFieldData::EntData &col_data) {
1268
1269#ifndef NDEBUG
1270 if (OpBase::locMat.size2() % SPACE_DIM)
1271 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1272 "Number of rows in matrix should be multiple of space dimensions");
1273#endif
1274
1275 // When we move to C++17 add if constexpr()
1276 if constexpr (COORDINATE_SYSTEM == POLAR || COORDINATE_SYSTEM == SPHERICAL)
1277 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1278 "%s coordiante not implemented",
1279 CoordinateTypesNames[COORDINATE_SYSTEM]);
1280
1281 auto t_w = this->getFTensor0IntegrationWeight();
1282 auto t_coords = this->getFTensor1CoordsAtGaussPts();
1283 size_t nb_base_functions_row = row_data.getN().size2();
1284 auto t_row_base = row_data.getFTensor0N();
1285 const double vol = this->getMeasure();
1286 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1287
1288 const double alpha =
1289 alphaConstant(t_coords(0), t_coords(1), t_coords(2)) * t_w * vol;
1290
1291 size_t rr = 0;
1292 auto t_m = getFTensor1FromPtr<SPACE_DIM>(OpBase::locMat.data().data());
1293
1294 // When we move to C++17 add if constexpr()
1295 if constexpr (COORDINATE_SYSTEM == CARTESIAN) {
1296 for (; rr != OpBase::nbRows; ++rr) {
1297 const double r_val = alpha * t_row_base;
1298 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1299 for (size_t cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1300 t_m(i) += r_val * t_col_diff_base(i);
1301 ++t_col_diff_base;
1302 ++t_m;
1303 }
1304 ++t_row_base;
1305 }
1306 }
1307
1308 // When we move to C++17 add if constexpr()
1309 if constexpr (COORDINATE_SYSTEM == CYLINDRICAL) {
1310 for (; rr != OpBase::nbRows; ++rr) {
1311 const double r_val = alpha * t_row_base;
1312 auto t_col_base = col_data.getFTensor0N(gg, 0);
1313 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1314 for (size_t cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1315 t_m(i) += r_val * t_col_diff_base(i);
1316 t_m(0) += (r_val / t_coords(0)) * t_col_base;
1317 ++t_col_base;
1318 ++t_col_diff_base;
1319 ++t_m;
1320 }
1321 ++t_row_base;
1322 }
1323 }
1324
1325 for (; rr < nb_base_functions_row; ++rr)
1326 ++t_row_base;
1327
1328 ++t_w;
1329 ++t_coords;
1330 }
1331
1333}
1334
1335template <int SPACE_DIM, typename OpBase>
1339 EntitiesFieldData::EntData &col_data) {
1341
1342 auto t_w = this->getFTensor0IntegrationWeight();
1343
1344 size_t nb_base_functions = row_data.getN().size2() / 3;
1345 auto t_row_base = row_data.getFTensor1N<3>();
1346 auto &mat = this->locMat;
1347 const double alpha_constant = alphaConstant();
1348 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1349
1350 const double alpha = alpha_constant * this->getMeasure() * t_w;
1351
1352 size_t rr = 0;
1353 for (; rr != OpBase::nbRows; ++rr) {
1354 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1355 for (size_t cc = 0; cc != OpBase::nbCols; ++cc) {
1356 mat(rr, cc) += alpha * t_row_base(i) * t_col_diff_base(i);
1357 ++t_col_diff_base;
1358 }
1359 ++t_row_base;
1360 }
1361 for (; rr < nb_base_functions; ++rr)
1362 ++t_row_base;
1363
1364 ++t_w;
1365 }
1366
1368}
1369
1370template <int SPACE_DIM, typename OpBase>
1374 EntitiesFieldData::EntData &col_data) {
1376
1377 auto t_w = this->getFTensor0IntegrationWeight();
1378
1379 size_t nb_base_functions = row_data.getN().size2();
1380 auto t_row_base = row_data.getFTensor0N();
1381
1382 auto get_t_vec = [&](const int rr) {
1383 std::array<double *, SPACE_DIM> ptrs;
1384 for (auto i = 0; i != SPACE_DIM; ++i)
1385 ptrs[i] = &OpBase::locMat(rr + i, 0);
1387 };
1388
1389 const double alpha_constant = alphaConstant();
1390 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1391
1392 const double alpha = alpha_constant * this->getMeasure() * t_w;
1393
1394 size_t rr = 0;
1395 for (; rr != OpBase::nbRows / SPACE_DIM; ++rr) {
1396 auto t_vec = get_t_vec(SPACE_DIM * rr);
1397 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1398 for (size_t cc = 0; cc != OpBase::nbCols; ++cc) {
1399 t_vec(i) += alpha * t_row_base * t_col_diff_base(i);
1400 ++t_col_diff_base;
1401 ++t_vec;
1402 }
1403 ++t_row_base;
1404 }
1405 for (; rr < nb_base_functions; ++rr)
1406 ++t_row_base;
1407
1408 ++t_w;
1409 }
1410
1412}
1413
1414template <int SPACE_DIM, typename OpBase>
1417 EntitiesFieldData::EntData &col_data) {
1419
1420 auto t_w = this->getFTensor0IntegrationWeight();
1421 auto t_coords = this->getFTensor1CoordsAtGaussPts();
1422
1423 size_t nb_base_functions = row_data.getN().size2() / 3;
1424 auto t_row_base = row_data.getFTensor1N<3>();
1425 const double alpha_constant = alphaConstant() * this->getMeasure();
1426 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1427
1428 const double alpha =
1429 alpha_constant * betaCoeff(t_coords(0), t_coords(1), t_coords(2)) * t_w;
1430
1431 size_t rr = 0;
1432 for (; rr != OpBase::nbRows / SPACE_DIM; ++rr) {
1433 auto t_mat_diag = getFTensor1FromArrayDiag<SPACE_DIM, SPACE_DIM>(
1434 this->locMat, SPACE_DIM * rr);
1435 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1436
1437 for (size_t cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1438 t_mat_diag(i) += alpha * t_row_base(j) * t_col_diff_base(j);
1439 ++t_col_diff_base;
1440 ++t_mat_diag;
1441 }
1442
1443 ++t_row_base;
1444 }
1445 for (; rr < nb_base_functions; ++rr)
1446 ++t_row_base;
1447
1448 ++t_w;
1449 ++t_coords;
1450 }
1451
1453}
1454
1455} // namespace MoFEM
1456
1457#endif //__BILINEAR_FORMS_INTEGRATORS_HPP__
constexpr double a
constexpr int SPACE_DIM
constexpr int FIELD_DIM
static const char *const CoordinateTypesNames[]
Coordinate system names.
Definition: definitions.h:108
#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
constexpr int BASE_DIM
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:308
FTensor::Index< 'i', SPACE_DIM > i
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: Common.hpp:10
double scalar_fun_one(const double, const double, const double)
boost::function< double()> ConstantFun
Constant function type.
constexpr IntegrationType I
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.
const VectorFieldEntities & getFieldEntities() const
get field entities
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
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)
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)
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, ScalarFun beta=scalar_fun_one)
OpGradTensorGradImpl(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< MatrixDouble > mat_D, ScalarFun beta=scalar_fun_one)
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)
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)
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)
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)
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, const bool only_transpose=false)
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)
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)
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)
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)
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)