v0.13.2
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 BASE_DIM, typename OpBase>
94struct OpMassImpl<3, BASE_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, 9, 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 <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S, IntegrationType I,
129 typename OpBase>
131
132template <int SPACE_DIM, int S, typename OpBase>
134 : public OpBase {
136 OpGradSymTensorGradImpl(const std::string row_field_name,
137 const std::string col_field_name,
138 boost::shared_ptr<MatrixDouble> mat_D,
139 boost::shared_ptr<Range> ents_ptr = nullptr,
141 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
142 matD(mat_D), betaCoeff(beta) {
143 if (row_field_name == col_field_name)
144 this->sYmm = true;
145 }
146
147protected:
148 boost::shared_ptr<MatrixDouble> matD;
150 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
152};
153
154template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S, IntegrationType I,
155 typename OpBase>
157
158template <int SPACE_DIM, int S, typename OpBase>
160 : public OpBase {
162 OpGradTensorGradImpl(const std::string row_field_name,
163 const std::string col_field_name,
164 boost::shared_ptr<MatrixDouble> mat_D,
166 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL), matD(mat_D),
167 betaCoeff(beta) {}
168
169protected:
170 boost::shared_ptr<MatrixDouble> matD;
172 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
174};
175
176template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S, IntegrationType I,
177 typename OpBase>
179
180template <int SPACE_DIM, int S, typename OpBase>
182 : public OpBase {
184 OpGradGradSymTensorGradGradImpl(const std::string row_field_name,
185 const std::string col_field_name,
186 boost::shared_ptr<MatrixDouble> mat_D,
187 boost::shared_ptr<Range> ents_ptr = nullptr)
188 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
189 matD(mat_D) {
190 if (row_field_name == col_field_name)
191 this->sYmm = true;
192 }
193
194protected:
195 boost::shared_ptr<MatrixDouble> matD;
196 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
198};
199
200template <int SPACE_DIM, IntegrationType I, typename OpBase>
202
203template <int SPACE_DIM, typename OpBase>
205 OpMixDivTimesScalarImpl(const std::string row_field_name,
206 const std::string col_field_name,
207 ConstantFun alpha_fun,
208 const bool assemble_transpose = false,
209 const bool only_transpose = false)
210 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
211 alphaConstant(alpha_fun) {
212 this->assembleTranspose = assemble_transpose;
213 this->onlyTranspose = only_transpose;
214 }
215
216protected:
219 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
221};
222
223template <int SPACE_DIM, IntegrationType I, typename OpBase>
225
226template <int SPACE_DIM, typename OpBase>
228 OpMixDivTimesVecImpl(const std::string row_field_name,
229 const std::string col_field_name, ConstantFun alpha_fun,
230 const bool assemble_transpose = false,
231 const bool only_transpose = false)
232 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
233 alphaConstant(alpha_fun) {
234 this->assembleTranspose = assemble_transpose;
235 this->onlyTranspose = only_transpose;
236 }
237
238protected:
241 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
243};
244
245template <int SPACE_DIM, IntegrationType I, typename OpBase,
246 CoordinateTypes COORDINATE_SYSTEM>
248
249template <int SPACE_DIM, typename OpBase, CoordinateTypes COORDINATE_SYSTEM>
250struct OpMixScalarTimesDivImpl<SPACE_DIM, GAUSS, OpBase, COORDINATE_SYSTEM>
251 : public OpBase {
253 const std::string row_field_name, const std::string col_field_name,
254 ScalarFun alpha_fun = [](double, double, double) constexpr { return 1; },
255 const bool assemble_transpose = false, const bool only_transpose = false)
256 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
257 alphaConstant(alpha_fun) {
258 this->assembleTranspose = assemble_transpose;
259 this->onlyTranspose = only_transpose;
260 this->sYmm = false;
261 }
262
263protected:
266 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
268};
269
270template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, IntegrationType I,
271 typename OpBase>
273
274template <int SPACE_DIM, typename OpBase>
276 : public OpBase {
277 OpMixVectorTimesGradImpl(const std::string row_field_name,
278 const std::string col_field_name,
279 ConstantFun alpha_fun,
280 const bool assemble_transpose = false,
281 const bool only_transpose = false)
282 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
283 alphaConstant(alpha_fun) {
284 this->assembleTranspose = assemble_transpose;
285 this->onlyTranspose = only_transpose;
286 }
287
288protected:
291 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
293};
294
295template <int SPACE_DIM, typename OpBase>
297 : public OpBase {
298 OpMixVectorTimesGradImpl(const std::string row_field_name,
299 const std::string col_field_name,
300 ConstantFun alpha_fun,
301 const bool assemble_transpose = false,
302 const bool only_transpose = false)
303 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
304 alphaConstant(alpha_fun) {
305 this->assembleTranspose = assemble_transpose;
306 this->onlyTranspose = only_transpose;
307 }
308
309protected:
312 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
314};
315
316template <int SPACE_DIM, IntegrationType I, typename OpBase>
318
319template <int SPACE_DIM, typename OpBase>
321 OpMixTensorTimesGradImpl(const std::string row_field_name,
322 const std::string col_field_name,
323 ConstantFun alpha_fun,
324 const bool assemble_transpose = false,
325 const bool only_transpose = false)
326 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
327 alphaConstant(alpha_fun) {
328 this->assembleTranspose = assemble_transpose;
329 this->onlyTranspose = only_transpose;
330 }
331
332protected:
336 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
338};
339
340template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, IntegrationType I,
341 typename OpBase>
343
344template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, IntegrationType I,
345 typename OpBase>
347
348template <int SPACE_DIM, typename OpBase>
350 : public OpBase {
352 const std::string field_name_row, const std::string field_name_col,
353 boost::shared_ptr<MatrixDouble> y_grad_ptr,
354 ConstantFun alpha_fun = []() { return 1; })
355 : OpBase(field_name_row, field_name_col, OpBase::OPROWCOL),
356 yGradPtr(y_grad_ptr), alphaConstant(alpha_fun) {
357
358 this->assembleTranspose = false;
359 this->onlyTranspose = false;
360 this->sYmm = false;
361 }
362
363protected:
364 boost::shared_ptr<MatrixDouble> yGradPtr;
366 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
368};
369
370template <int SPACE_DIM, typename OpBase>
372 : public OpBase {
374 const std::string field_name_row, const std::string field_name_col,
375 boost::shared_ptr<MatrixDouble> u_ptr,
376 ConstantFun alpha_fun = []() { return 1; })
377 : OpBase(field_name_row, field_name_col, OpBase::OPROWCOL), uPtr(u_ptr),
378 alphaConstant(alpha_fun) {
379
380 this->assembleTranspose = false;
381 this->onlyTranspose = false;
382 this->sYmm = false;
383 }
384
385protected:
387 boost::shared_ptr<MatrixDouble> uPtr;
388 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
390};
391
392template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
394 : public OpBase {
396 const std::string field_name_row, const std::string field_name_col,
397 boost::shared_ptr<MatrixDouble> y_grad_ptr,
398 ConstantFun alpha_fun = []() { return 1; })
399 : OpBase(field_name_row, field_name_col, OpBase::OPROWCOL),
400 yGradPtr(y_grad_ptr), alphaConstant(alpha_fun) {
401
402 this->assembleTranspose = false;
403 this->onlyTranspose = false;
404 this->sYmm = false;
405 }
406
407protected:
409 boost::shared_ptr<MatrixDouble> yGradPtr;
410 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
412};
413
414template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
416 : public OpBase {
418 const std::string field_name_row, const std::string field_name_col,
419 boost::shared_ptr<MatrixDouble> u_ptr,
420 ConstantFun alpha_fun = []() { return 1; })
421 : OpBase(field_name_row, field_name_col, OpBase::OPROWCOL), uPtr(u_ptr),
422 alphaConstant(alpha_fun) {
423
424 this->assembleTranspose = false;
425 this->onlyTranspose = false;
426 this->sYmm = false;
427 }
428
429protected:
431 boost::shared_ptr<MatrixDouble> uPtr;
432 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
434};
435
436/**
437 * @brief Bilinear integrator form
438 * @ingroup mofem_forms
439 *
440 * @tparam EleOp
441 * @tparam A
442 * @tparam I
443 */
444template <typename EleOp>
445template <AssemblyType A>
446template <IntegrationType I>
448
449 /**
450 * @brief Integrate \f$(v_{,i},\beta(\mathbf{x}) u_{,j}))_\Omega\f$
451 * @ingroup mofem_forms
452 *
453 * @tparam SPACE_DIM
454 */
455 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM>
457
458 /**
459 * @brief Integrate \f$(v_i,\beta(\mathbf{x}) u_j)_\Omega\f$
460 * @ingroup mofem_forms
461 *
462 * @tparam
463 */
464 template <int BASE_DIM, int FIELD_DIM>
466
467 /**
468 * @brief Integrate \f$(v_k,D_{ijkl} u_{,l})_\Omega\f$
469 *
470 * \note \f$D_{ijkl}\f$ is * tensor with minor & major symmetry
471 *
472 * @ingroup mofem_forms
473 *
474 * @tparam SPACE_DIM
475 * @tparam S
476 */
477 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S = 0>
480
481 /**
482 * @brief Integrate \f$(v_{,ij},D_{ijkl} u_{,kl})_\Omega\f$
483 *
484 * \note \f$D_{ijkl}\f$ is * tensor with no symmetries
485 *
486 * @ingroup mofem_forms
487 *
488 * @tparam SPACE_DIM
489 * @tparam S
490 */
491 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S = 0>
494 OpBase>;
495
496 /**
497 * @brief Integrate \f$(v_{,j},D_{ijkl} u_{,l})_\Omega\f$
498 *
499 * \note \f$D_{ijkl}\f$ is * tensor with no symmetries
500 *
501 * @ingroup mofem_forms
502 *
503 * @tparam SPACE_DIM
504 * @tparam S
505 */
506 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S = 0>
509
510 /**
511 * @brief Integrate \f$(\lambda_{i,i},u)_\Omega\f$
512 *
513 * @tparam SPACE_DIM
514 */
515 template <int SPACE_DIM>
517
518 /**
519 * @brief Integrate \f$(\lambda_{ij,j},u_{i})_\Omega\f$
520 *
521 * @tparam SPACE_DIM
522 */
523 template <int SPACE_DIM>
525
526 /**
527 * @brief Integrate \f$(\lambda,u_{i,i})_\Omega\f$
528 *
529 * @tparam SPACE_DIM
530 */
531 template <int SPACE_DIM, CoordinateTypes COORDINATE_SYSTEM = CARTESIAN>
534
535 /**
536 * @brief Integrate \f$(\lambda_{i},u_{,j})_\Omega\f$
537 *
538 * @tparam SPACE_DIM
539 */
540 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM>
543
544 /**
545 * @brief Integrate \f$(\lambda_{ij},u_{i,j})_\Omega\f$
546 *
547 * @tparam SPACE_DIM
548 */
549 template <int SPACE_DIM>
551
552 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM>
555
556 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM>
559};
560
561template <int SPACE_DIM, typename OpBase>
564 EntitiesFieldData::EntData &col_data) {
566
567 // get element volume
568 const double vol = OpBase::getMeasure();
569 // get integration weights
570 auto t_w = OpBase::getFTensor0IntegrationWeight();
571 // get base function gradient on rows
572 auto t_row_grad = row_data.getFTensor1DiffN<SPACE_DIM>();
573 // get coordinate at integration points
574 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
575
576 // loop over integration points
577 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
578 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
579 // take into account Jacobian
580 const double alpha = t_w * beta;
581 // loop over ros base functions
582 int rr = 0;
583 for (; rr != OpBase::nbRows; rr++) {
584 // get column base functions gradient at gauss point gg
585 auto t_col_grad = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
586 // loop over columns
587 for (int cc = 0; cc != OpBase::nbCols; cc++) {
588 // calculate element of local matrix
589 OpBase::locMat(rr, cc) += alpha * (t_row_grad(i) * t_col_grad(i));
590 ++t_col_grad; // move to another gradient of base function
591 // on column
592 }
593 ++t_row_grad; // move to another element of gradient of base
594 // function on row
595 }
596 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
597 ++t_row_grad;
598
599 ++t_coords;
600 ++t_w; // move to another integration weight
601 }
603}
604
605template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
609 EntitiesFieldData::EntData &col_data) {
611
612 // get element volume
613 const double vol = OpBase::getMeasure();
614 // get integration weights
615 auto t_w = OpBase::getFTensor0IntegrationWeight();
616 // get base function gradient on rows
617 auto t_row_grad = row_data.getFTensor1DiffN<SPACE_DIM>();
618 // get coordinate at integration points
619 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
620
623
624 auto get_t_vec = [&](const int rr) {
625 std::array<double *, FIELD_DIM> ptrs;
626 for (auto i = 0; i != FIELD_DIM; ++i)
627 ptrs[i] = &OpBase::locMat(rr + i, i);
629 ptrs);
630 };
631
632 // loop over integration points
633 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
634 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
635 // take into account Jacobian
636 const double alpha = t_w * beta;
637 // loop over ros base functions
638 int rr = 0;
639 for (; rr != OpBase::nbRows / FIELD_DIM; rr++) {
640 // get diag vec
641 auto t_vec = get_t_vec(rr * FIELD_DIM);
642 // get column base functions gradient at gauss point gg
643 auto t_col_grad = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
644 // loop over columns
645 for (int cc = 0; cc != OpBase::nbCols / FIELD_DIM; cc++) {
646 // calculate element of local matrix
647 t_vec(j) += alpha * (t_row_grad(i) * t_col_grad(i));
648 ++t_col_grad; // move to another gradient of base function
649 // on column
650 ++t_vec;
651 }
652 ++t_row_grad; // move to another element of gradient of base
653 // function on row
654 }
655 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
656 ++t_row_grad;
657
658 ++t_coords;
659 ++t_w; // move to another integration weight
660 }
662}
663
664template <typename OpBase>
667 EntitiesFieldData::EntData &col_data) {
669
670#ifndef NDEBUG
671 auto log_error = [&]() {
672 MOFEM_LOG("SELF", Sev::error) << "Row side " << OpBase::rowSide << " "
673 << CN::EntityTypeName(OpBase::rowType);
674 MOFEM_LOG("SELF", Sev::error) << "Col side " << OpBase::colSide << " "
675 << CN::EntityTypeName(OpBase::colType);
676 };
677
678 if (row_data.getN().size2() < OpBase::nbRows) {
679 log_error();
680 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
681 "Wrong number of base functions on rows %d < %d",
682 row_data.getN().size2(), OpBase::nbRows);
683 }
684 if (col_data.getN().size2() < OpBase::nbCols) {
685 log_error();
686 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
687 "Wrong number of base functions on cols %d < %d",
688 col_data.getN().size2(), OpBase::nbCols);
689 }
690 if (row_data.getN().size1() != OpBase::nbIntegrationPts) {
691 log_error();
692 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
693 "Wrong number of integration points on rows %d != %d",
694 row_data.getN().size1(), OpBase::nbIntegrationPts);
695 }
696 if (col_data.getN().size1() != OpBase::nbIntegrationPts) {
697 log_error();
698 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
699 "Wrong number of integration points on cols %d != %d",
700 col_data.getN().size1(), OpBase::nbIntegrationPts);
701 }
702#endif
703
704 // get element volume
705 const double vol = OpBase::getMeasure();
706 // get integration weights
707 auto t_w = OpBase::getFTensor0IntegrationWeight();
708 // get base function gradient on rows
709 auto t_row_base = row_data.getFTensor0N();
710 // get coordinate at integration points
711 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
712 // loop over integration points
713 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
714 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
715 // take into account Jacobian
716 const double alpha = t_w * beta;
717 // loop over rows base functions
718 auto a_mat_ptr = &*OpBase::locMat.data().begin();
719 int rr = 0;
720 for (; rr != OpBase::nbRows; rr++) {
721 // get column base functions gradient at gauss point gg
722 auto t_col_base = col_data.getFTensor0N(gg, 0);
723 // loop over columns
724 for (int cc = 0; cc != OpBase::nbCols; cc++) {
725 // calculate element of local matrix
726 *a_mat_ptr += alpha * (t_row_base * t_col_base);
727 ++t_col_base;
728 ++a_mat_ptr;
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) {
745 // get element volume
746 const double vol = OpBase::getMeasure();
747 // get integration weights
748 auto t_w = OpBase::getFTensor0IntegrationWeight();
749 // get base function gradient on rows
750 auto t_row_base = row_data.getFTensor0N();
751 // get coordinate at integration points
752 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
753
755 auto get_t_vec = [&](const int rr) {
756 std::array<double *, FIELD_DIM> ptrs;
757 for (auto i = 0; i != FIELD_DIM; ++i)
758 ptrs[i] = &OpBase::locMat(rr + i, i);
760 ptrs);
761 };
762
763 // loop over integration points
764 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
765 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
766 // take into account Jacobian
767 const double alpha = t_w * beta;
768 // loop over rows base functions
769 int rr = 0;
770 for (; rr != OpBase::nbRows / FIELD_DIM; rr++) {
771 // get column base functions gradient at gauss point gg
772 auto t_col_base = col_data.getFTensor0N(gg, 0);
773 // get mat vec
774 auto t_vec = get_t_vec(FIELD_DIM * rr);
775 // loop over columns
776 for (int cc = 0; cc != OpBase::nbCols / FIELD_DIM; cc++) {
777 // calculate element of local matrix
778 t_vec(i) += alpha * (t_row_base * t_col_base);
779 ++t_col_base;
780 ++t_vec;
781 }
782 ++t_row_base;
783 }
784 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
785 ++t_row_base;
786 ++t_coords;
787 ++t_w; // move to another integration weight
788 }
790};
791
792template <int BASE_DIM, typename OpBase>
795 EntitiesFieldData::EntData &col_data) {
798 size_t nb_base_functions = row_data.getN().size2() / 3;
799 // // get element volume
800 const double vol = OpBase::getMeasure();
801 // get integration weights
802 auto t_w = OpBase::getFTensor0IntegrationWeight();
803 // get base function gradient on rows
804 auto t_row_base = row_data.getFTensor1N<3>();
805 // get coordinate at integration points
806 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
807 // loop over integration points
808 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
809 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
810 // take into account Jacobian
811 const double alpha = t_w * beta;
812 // loop over rows base functions
813 auto a_mat_ptr = &*OpBase::locMat.data().begin();
814 int rr = 0;
815 for (; rr != OpBase::nbRows; rr++) {
816 // get column base functions gradient at gauss point gg
817 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
818 // loop over columns
819 for (int cc = 0; cc != OpBase::nbCols; cc++) {
820 // calculate element of local matrix
821 (*a_mat_ptr) += alpha * (t_row_base(i) * t_col_base(i));
822 ++t_col_base;
823 ++a_mat_ptr;
824 }
825 ++t_row_base;
826 }
827 for (; rr < nb_base_functions; ++rr)
828 ++t_row_base;
829 ++t_coords;
830 ++t_w; // move to another integration weight
831 }
833};
834
835template <typename OpBase>
838 EntitiesFieldData::EntData &col_data) {
842 auto get_t_vec = [&](const int rr) {
843 std::array<double *, 3> ptrs;
844 for (auto i = 0; i != 3; ++i)
845 ptrs[i] = &OpBase::locMat(rr + i, i);
847 };
848 size_t nb_base_functions = row_data.getN().size2() / 3;
849 // // get element volume
850 const double vol = OpBase::getMeasure();
851 // get integration weights
852 auto t_w = OpBase::getFTensor0IntegrationWeight();
853 // get base function gradient on rows
854 auto t_row_base = row_data.getFTensor1N<3>();
855 // get coordinate at integration points
856 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
857 // loop over integration points
858 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
859 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
860 // take into account Jacobian
861 const double alpha = t_w * beta;
862 // loop over rows base functions
863 int rr = 0;
864 for (; rr != OpBase::nbRows / 3; rr++) {
865 // get column base functions gradient at gauss point gg
866 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
867 auto t_vec = get_t_vec(3 * rr);
868 // loop over columns
869 for (int cc = 0; cc != OpBase::nbCols / 3; cc++) {
870 // calculate element of local matrix
871 t_vec(i) += alpha * (t_row_base(k) * t_col_base(k));
872 ++t_col_base;
873 ++t_vec;
874 }
875 ++t_row_base;
876 }
877 for (; rr < nb_base_functions; ++rr)
878 ++t_row_base;
879 ++t_coords;
880 ++t_w; // move to another integration weight
881 }
883};
884
885template <int SPACE_DIM, int S, typename OpBase>
889 EntitiesFieldData::EntData &col_data) {
891
892 const size_t nb_row_dofs = row_data.getIndices().size();
893 const size_t nb_col_dofs = col_data.getIndices().size();
894
895 if (nb_row_dofs && nb_col_dofs) {
896
897 const bool diag = (row_data.getFieldEntities()[0]->getLocalUniqueId() ==
898 col_data.getFieldEntities()[0]->getLocalUniqueId());
899
904
905 // get element volume
906 double vol = OpBase::getMeasure();
907
908 // get intergrayion weights
909 auto t_w = OpBase::getFTensor0IntegrationWeight();
910
911 // get derivatives of base functions on rows
912 auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
913
914 // Elastic stiffness tensor (4th rank tensor with minor and major
915 // symmetry)
916 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, S>(*matD);
917
918 // get coordinate at integration points
919 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
920
921 // iterate over integration points
922 for (int gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
923
924 // calculate scalar weight times element volume
925 double a = t_w * vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
926
927 // iterate over row base functions
928 int rr = 0;
929 for (; rr != OpBase::nbRows / SPACE_DIM; ++rr) {
930
931 // get sub matrix for the row
932 auto t_m = OpBase::template getLocMat<SPACE_DIM>(SPACE_DIM * rr);
933
935 // I mix up the indices here so that it behaves like a
936 // Dg. That way I don't have to have a separate wrapper
937 // class Christof_Expr, which simplifies things.
938 t_rowD(l, j, k) = t_D(i, j, k, l) * (a * t_row_diff_base(i));
939
940 // get derivatives of base functions for columns
941 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
942
943 // iterate column base functions
944 const auto nb_cols = (diag) ? rr : OpBase::nbCols / SPACE_DIM - 1;
945 for (int cc = 0; cc <= nb_cols; ++cc) {
946
947 // integrate block local stiffens matrix
948 t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base(k);
949
950 // move to next column base function
951 ++t_col_diff_base;
952
953 // move to next block of local stiffens matrix
954 ++t_m;
955 }
956
957 // move to next row base function
958 ++t_row_diff_base;
959 }
960
961 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
962 ++t_row_diff_base;
963
964 // move to next integration weight
965 ++t_w;
966 ++t_D;
967 ++t_coords;
968 }
969
970 // Copy symmetry
971 if (diag) {
972 for (int rr = 0; rr != OpBase::nbRows / SPACE_DIM - 1; ++rr) {
973 auto t_m_rr = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
974 this->locMat, SPACE_DIM * rr, SPACE_DIM * (rr + 1));
975 auto t_m_cc = getFTensor2FromArray<SPACE_DIM, SPACE_DIM>(
976 this->locMat, SPACE_DIM * (rr + 1), SPACE_DIM * rr,
978 for (int cc = rr + 1; cc < OpBase::nbCols / SPACE_DIM; ++cc) {
979 t_m_rr(i, j) = t_m_cc(j, i);
980 ++t_m_rr;
981 ++t_m_cc;
982 }
983 }
984 }
985
986 }
987
989}
990
991template <int SPACE_DIM, int S, typename OpBase>
995 EntitiesFieldData::EntData &col_data) {
997
999
1004
1005 auto &row_hessian = row_data.getN(BaseDerivatives::SecondDerivative);
1006 auto &col_hessian = col_data.getN(BaseDerivatives::SecondDerivative);
1007
1008#ifndef NDEBUG
1009 if (row_hessian.size1() != OpBase::nbIntegrationPts) {
1010 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1011 "Wrong number of integration pts (%d != %d)",
1012 row_hessian.size1(), OpBase::nbIntegrationPts);
1013 }
1014 if (row_hessian.size2() !=
1016 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1017 "Wrong number of base functions (%d != %d)",
1018 row_hessian.size2() / (SPACE_DIM * SPACE_DIM),
1020 }
1021 if (row_hessian.size2() < OpBase::nbRows * SPACE_DIM * SPACE_DIM) {
1022 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1023 "Wrong number of base functions (%d < %d)", row_hessian.size2(),
1025 }
1026 if (col_hessian.size1() != OpBase::nbIntegrationPts) {
1027 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1028 "Wrong number of integration pts (%d != %d)",
1029 col_hessian.size1(), OpBase::nbIntegrationPts);
1030 }
1031 if (col_hessian.size2() < OpBase::nbCols * SPACE_DIM * SPACE_DIM) {
1032 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1033 "Wrong number of base functions (%d < %d)", col_hessian.size2(),
1035 }
1036#endif
1037
1038 // get element volume
1039 double vol = OpBase::getMeasure();
1040
1041 // get intergrayion weights
1042 auto t_w = OpBase::getFTensor0IntegrationWeight();
1043
1044 auto t_row_diff2 = getFTensor2SymmetricLowerFromPtr<SPACE_DIM>(
1045 &*row_hessian.data().begin());
1046
1047 // Elastic stiffness tensor (4th rank tensor with minor and major
1048 // symmetry)
1049 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, S>(*matD);
1050
1051 // iterate over integration points
1052 for (int gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1053
1054 // calculate scalar weight times element volume
1055 double a = t_w * vol;
1056
1057 // get sub matrix for the row
1058 auto m_ptr = &*OpBase::locMat.data().begin();
1059
1060 // iterate over row base functions
1061 int rr = 0;
1062 for (; rr != OpBase::nbRows; ++rr) {
1063
1065 t_rowD(k, l) = t_D(i, j, k, l) * (a * t_row_diff2(i, j));
1066
1067 // get derivatives of base functions for columns
1068 auto t_col_diff2 =
1069 getFTensor2SymmetricLowerFromPtr<SPACE_DIM>(&col_hessian(gg, 0));
1070
1071 // iterate column base functions
1072 for (int cc = 0; cc != OpBase::nbCols; ++cc) {
1073
1074 // integrate block local stiffens matrix
1075 *m_ptr += t_rowD(i, j) * t_col_diff2(i, j);
1076
1077 // move to next column base function
1078 ++t_col_diff2;
1079
1080 // move to next block of local stiffens matrix
1081 ++m_ptr;
1082 }
1083
1084 // move to next row base function
1085 ++t_row_diff2;
1086 }
1087
1088 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
1089 ++t_row_diff2;
1090
1091 // move to next integration weight
1092 ++t_w;
1093 ++t_D;
1094 }
1095 }
1096
1098}
1099
1100template <int SPACE_DIM, int S, typename OpBase>
1104 EntitiesFieldData::EntData &col_data) {
1106
1107 const size_t nb_row_dofs = row_data.getIndices().size();
1108 const size_t nb_col_dofs = col_data.getIndices().size();
1109
1110 if (nb_row_dofs && nb_col_dofs) {
1111
1116
1117 // get element volume
1118 double vol = OpBase::getMeasure();
1119
1120 // get intergrayion weights
1121 auto t_w = OpBase::getFTensor0IntegrationWeight();
1122
1123 // get derivatives of base functions on rows
1124 auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
1125
1126 // stiffness tensor (4th rank tensor)
1127 auto t_D =
1128 getFTensor4FromMat<SPACE_DIM, SPACE_DIM, SPACE_DIM, SPACE_DIM, S>(
1129 *matD);
1130
1131 // get coordinate at integration points
1132 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
1133
1134 // iterate over integration points
1135 for (int gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1136
1137 // calculate scalar weight times element volume
1138 double a = t_w * vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
1139
1140 // iterate over row base functions
1141 int rr = 0;
1142 for (; rr != OpBase::nbRows / SPACE_DIM; ++rr) {
1143
1144 // get sub matrix for the row
1145 auto t_m = OpBase::template getLocMat<SPACE_DIM>(SPACE_DIM * rr);
1146
1147 // get derivatives of base functions for columns
1148 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1149
1150 // iterate column base functions
1151 for (int cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1152
1153 // integrate block local stiffens matrix
1154 t_m(i, k) +=
1155 a * (t_D(i, j, k, l) * (t_row_diff_base(j) * t_col_diff_base(l)));
1156
1157 // move to next column base function
1158 ++t_col_diff_base;
1159
1160 // move to next block of local stiffens matrix
1161 ++t_m;
1162 }
1163
1164 // move to next row base function
1165 ++t_row_diff_base;
1166 }
1167
1168 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
1169 ++t_row_diff_base;
1170
1171 // move to next integration weight
1172 ++t_w;
1173 ++t_D;
1174 ++t_coords;
1175 }
1176 }
1177
1179}
1180
1181template <int SPACE_DIM, typename OpBase>
1184 EntitiesFieldData::EntData &col_data) {
1186
1187 auto t_w = this->getFTensor0IntegrationWeight();
1188
1189 size_t nb_base_functions = row_data.getN().size2() / 3;
1190 auto t_row_diff_base = row_data.getFTensor2DiffN<3, SPACE_DIM>();
1191
1192 const double alpha_constant = alphaConstant();
1193
1194 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1195 const double alpha = alpha_constant * this->getMeasure() * t_w;
1196
1197 size_t rr = 0;
1198 for (; rr != OpBase::nbRows; ++rr) {
1199 const double t_row_div_base = t_row_diff_base(i, i);
1200 auto t_col_base = col_data.getFTensor0N(gg, 0);
1201 for (size_t cc = 0; cc != OpBase::nbCols; ++cc) {
1202 this->locMat(rr, cc) += alpha * t_row_div_base * t_col_base;
1203 ++t_col_base;
1204 }
1205 ++t_row_diff_base;
1206 }
1207 for (; rr < nb_base_functions; ++rr)
1208 ++t_row_diff_base;
1209
1210 ++t_w;
1211 }
1212
1214}
1215
1216template <int SPACE_DIM, typename OpBase>
1219 EntitiesFieldData::EntData &col_data) {
1221
1222 auto t_w = this->getFTensor0IntegrationWeight();
1223
1224 size_t nb_base_functions = row_data.getN().size2() / 3;
1225 auto t_row_diff_base = row_data.getFTensor2DiffN<3, SPACE_DIM>();
1226 const double alpha_constant = alphaConstant();
1227 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1228
1229 const double alpha = alpha_constant * this->getMeasure() * t_w;
1230
1231 size_t rr = 0;
1232 for (; rr != OpBase::nbRows / SPACE_DIM; ++rr) {
1233 auto t_mat_diag = getFTensor1FromArrayDiag<SPACE_DIM, SPACE_DIM>(
1234 this->locMat, SPACE_DIM * rr);
1235 const double t_row_div_base = t_row_diff_base(i, i);
1236 auto t_col_base = col_data.getFTensor0N(gg, 0);
1237
1238 for (size_t cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1239 t_mat_diag(i) += alpha * t_row_div_base * t_col_base;
1240 ++t_col_base;
1241 ++t_mat_diag;
1242 }
1243
1244 ++t_row_diff_base;
1245 }
1246 for (; rr < nb_base_functions; ++rr)
1247 ++t_row_diff_base;
1248
1249 ++t_w;
1250 }
1251
1253}
1254
1255template <int SPACE_DIM, typename OpBase, CoordinateTypes COORDINATE_SYSTEM>
1259 EntitiesFieldData::EntData &col_data) {
1261
1262#ifndef NDEBUG
1263 if (OpBase::locMat.size2() % SPACE_DIM)
1264 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1265 "Number of rows in matrix should be multiple of space dimensions");
1266#endif
1267
1268 // When we move to C++17 add if constexpr()
1269 if constexpr (COORDINATE_SYSTEM == POLAR || COORDINATE_SYSTEM == SPHERICAL)
1270 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1271 "%s coordiante not implemented",
1272 CoordinateTypesNames[COORDINATE_SYSTEM]);
1273
1274 auto t_w = this->getFTensor0IntegrationWeight();
1275 auto t_coords = this->getFTensor1CoordsAtGaussPts();
1276 size_t nb_base_functions_row = row_data.getN().size2();
1277 auto t_row_base = row_data.getFTensor0N();
1278 const double vol = this->getMeasure();
1279 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1280
1281 const double alpha =
1282 alphaConstant(t_coords(0), t_coords(1), t_coords(2)) * t_w * vol;
1283
1284 size_t rr = 0;
1285 auto t_m = getFTensor1FromPtr<SPACE_DIM>(OpBase::locMat.data().data());
1286
1287 // When we move to C++17 add if constexpr()
1288 if constexpr (COORDINATE_SYSTEM == CARTESIAN) {
1289 for (; rr != OpBase::nbRows; ++rr) {
1290 const double r_val = alpha * t_row_base;
1291 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1292 for (size_t cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1293 t_m(i) += r_val * t_col_diff_base(i);
1294 ++t_col_diff_base;
1295 ++t_m;
1296 }
1297 ++t_row_base;
1298 }
1299 }
1300
1301 // When we move to C++17 add if constexpr()
1302 if constexpr (COORDINATE_SYSTEM == CYLINDRICAL) {
1303 for (; rr != OpBase::nbRows; ++rr) {
1304 const double r_val = alpha * t_row_base;
1305 auto t_col_base = col_data.getFTensor0N(gg, 0);
1306 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1307 for (size_t cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1308 t_m(i) += r_val * t_col_diff_base(i);
1309 t_m(0) += (r_val / t_coords(0)) * t_col_base;
1310 ++t_col_base;
1311 ++t_col_diff_base;
1312 ++t_m;
1313 }
1314 ++t_row_base;
1315 }
1316 }
1317
1318 for (; rr < nb_base_functions_row; ++rr)
1319 ++t_row_base;
1320
1321 ++t_w;
1322 ++t_coords;
1323 }
1324
1326}
1327
1328template <int SPACE_DIM, typename OpBase>
1332 EntitiesFieldData::EntData &col_data) {
1334
1335 auto t_w = this->getFTensor0IntegrationWeight();
1336
1337 size_t nb_base_functions = row_data.getN().size2() / 3;
1338 auto t_row_base = row_data.getFTensor1N<3>();
1339 auto &mat = this->locMat;
1340 const double alpha_constant = alphaConstant();
1341 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1342
1343 const double alpha = alpha_constant * this->getMeasure() * t_w;
1344
1345 size_t rr = 0;
1346 for (; rr != OpBase::nbRows; ++rr) {
1347 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1348 for (size_t cc = 0; cc != OpBase::nbCols; ++cc) {
1349 mat(rr, cc) += alpha * t_row_base(i) * t_col_diff_base(i);
1350 ++t_col_diff_base;
1351 }
1352 ++t_row_base;
1353 }
1354 for (; rr < nb_base_functions; ++rr)
1355 ++t_row_base;
1356
1357 ++t_w;
1358 }
1359
1361}
1362
1363template <int SPACE_DIM, typename OpBase>
1367 EntitiesFieldData::EntData &col_data) {
1369
1370 auto t_w = this->getFTensor0IntegrationWeight();
1371
1372 size_t nb_base_functions = row_data.getN().size2();
1373 auto t_row_base = row_data.getFTensor0N();
1374
1375 auto get_t_vec = [&](const int rr) {
1376 std::array<double *, SPACE_DIM> ptrs;
1377 for (auto i = 0; i != SPACE_DIM; ++i)
1378 ptrs[i] = &OpBase::locMat(rr + i, 0);
1380 };
1381
1382 const double alpha_constant = alphaConstant();
1383 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1384
1385 const double alpha = alpha_constant * this->getMeasure() * t_w;
1386
1387 size_t rr = 0;
1388 for (; rr != OpBase::nbRows / SPACE_DIM; ++rr) {
1389 auto t_vec = get_t_vec(SPACE_DIM * rr);
1390 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1391 for (size_t cc = 0; cc != OpBase::nbCols; ++cc) {
1392 t_vec(i) += alpha * t_row_base * t_col_diff_base(i);
1393 ++t_col_diff_base;
1394 ++t_vec;
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>
1410 EntitiesFieldData::EntData &col_data) {
1412
1413 auto t_w = this->getFTensor0IntegrationWeight();
1414
1415 size_t nb_base_functions = row_data.getN().size2() / 3;
1416 auto t_row_base = row_data.getFTensor1N<3>();
1417 const double alpha_constant = alphaConstant();
1418 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1419
1420 const double alpha = alpha_constant * this->getMeasure() * t_w;
1421
1422 size_t rr = 0;
1423 for (; rr != OpBase::nbRows / SPACE_DIM; ++rr) {
1424 auto t_mat_diag = getFTensor1FromArrayDiag<SPACE_DIM, SPACE_DIM>(
1425 this->locMat, SPACE_DIM * rr);
1426 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1427
1428 for (size_t cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1429 t_mat_diag(i) += alpha * t_row_base(j) * t_col_diff_base(j);
1430 ++t_col_diff_base;
1431 ++t_mat_diag;
1432 }
1433
1434 ++t_row_base;
1435 }
1436 for (; rr < nb_base_functions; ++rr)
1437 ++t_row_base;
1438
1439 ++t_w;
1440 }
1441
1443}
1444
1445template <int SPACE_DIM, typename OpBase>
1449 EntitiesFieldData::EntData &col_data) {
1451
1452 // get element volume
1453 const double vol = OpBase::getMeasure();
1454 // get integration weights
1455 auto t_w = OpBase::getFTensor0IntegrationWeight();
1456 // get base function gradient on rows
1457 auto t_row_base = row_data.getFTensor0N();
1458
1459 auto get_t_vec = [&](const int rr) {
1460 std::array<double *, SPACE_DIM> ptrs;
1461 for (auto i = 0; i != SPACE_DIM; ++i)
1462 ptrs[i] = &OpBase::locMat(rr, i);
1464 ptrs);
1465 };
1466
1467 auto t_grad_y = getFTensor1FromMat<SPACE_DIM>(*yGradPtr);
1469
1470 const double alpha_constant = alphaConstant();
1471 // loop over integration points
1472 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
1473 // take into account Jacobian
1474 const double alpha = t_w * vol * alpha_constant;
1475 // access local matrix
1476 auto t_vec = get_t_vec(0);
1477 // loop over rows base functions
1478 int rr = 0;
1479 for (; rr != OpBase::nbRows; rr++) {
1480 // get column base functions gradient at gauss point gg
1481 auto t_col_base = col_data.getFTensor0N(gg, 0);
1482 // loop over columns
1483 for (int cc = 0; cc != OpBase::nbCols / SPACE_DIM; cc++) {
1484 // calculate element of local matrix
1485 t_vec(i) += alpha * t_row_base * t_col_base * t_grad_y(i);
1486 ++t_col_base;
1487 ++t_vec;
1488 }
1489 ++t_row_base;
1490 }
1491 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
1492 ++t_row_base;
1493
1494 ++t_grad_y;
1495 ++t_w; // move to another integration weight
1496 }
1498};
1499
1500template <int SPACE_DIM, typename OpBase>
1504 EntitiesFieldData::EntData &col_data) {
1506
1507 // get element volume
1508 const double vol = OpBase::getMeasure();
1509 // get integration weights
1510 auto t_w = OpBase::getFTensor0IntegrationWeight();
1511 // get base function gradient on rows
1512 auto t_row_base = row_data.getFTensor0N();
1513
1514 auto t_u = getFTensor1FromMat<SPACE_DIM>(*uPtr);
1516 const double alpha_constant = alphaConstant();
1517 // loop over integration points
1518 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
1519 // take into account Jacobian
1520 const double alpha = t_w * vol * alpha_constant;
1521 // loop over rows base functions
1522 auto a_mat_ptr = &*OpBase::locMat.data().begin();
1523 int rr = 0;
1524 for (; rr != OpBase::nbRows; rr++) {
1525 // get column base functions gradient at gauss point gg
1526 auto t_diff_col_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1527 // loop over columns
1528 for (int cc = 0; cc != OpBase::nbCols; cc++) {
1529 // calculate element of local matrix
1530 (*a_mat_ptr) += alpha * t_row_base * t_diff_col_base(i) * t_u(i);
1531 ++t_diff_col_base;
1532 ++a_mat_ptr;
1533 }
1534 ++t_row_base;
1535 }
1536 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
1537 ++t_row_base;
1538
1539 ++t_u;
1540 ++t_w; // move to another integration weight
1541 }
1543};
1544
1545template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
1549 EntitiesFieldData::EntData &col_data) {
1551
1552 // get element volume
1553 const double vol = OpBase::getMeasure();
1554 // get integration weights
1555 auto t_w = OpBase::getFTensor0IntegrationWeight();
1556 // get base function gradient on rows
1557 auto t_row_base = row_data.getFTensor0N();
1558
1559 auto t_grad_y = getFTensor2FromMat<FIELD_DIM, SPACE_DIM>(*yGradPtr);
1562
1563 auto get_t_mat = [&](const int rr) {
1564 std::array<double *, FIELD_DIM * SPACE_DIM> ptrs;
1565 int s = 0;
1566 for (int j = 0; j != FIELD_DIM; ++j)
1567 for (auto i = 0; i != SPACE_DIM; ++i, ++s)
1568 ptrs[s] = &OpBase::locMat(rr + j, i);
1570 SPACE_DIM>(ptrs);
1571 };
1572
1573 const double alpha_constant = alphaConstant();
1574 // loop over integration points
1575 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
1576 // take into account Jacobian
1577 const double alpha = t_w * vol * alpha_constant;
1578
1579 // loop over rows base functions
1580 int rr = 0;
1581 for (; rr != OpBase::nbRows / FIELD_DIM; rr++) {
1582 // get column base functions gradient at gauss point gg
1583 auto t_col_base = col_data.getFTensor0N(gg, 0);
1584 // get mat
1585 auto t_mat = get_t_mat(FIELD_DIM * rr);
1586 // loop over columns
1587 for (int cc = 0; cc != OpBase::nbCols / SPACE_DIM; cc++) {
1588 // calculate element of local matrix
1589 t_mat(I, k) += (alpha * t_row_base * t_col_base) * t_grad_y(I, k);
1590 ++t_col_base;
1591 ++t_mat;
1592 }
1593 ++t_row_base;
1594 }
1595 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
1596 ++t_row_base;
1597
1598 ++t_grad_y;
1599 ++t_w; // move to another integration weight
1600 }
1602};
1603
1604template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
1608 EntitiesFieldData::EntData &col_data) {
1610
1611 // get element volume
1612 const double vol = OpBase::getMeasure();
1613 // get integration weights
1614 auto t_w = OpBase::getFTensor0IntegrationWeight();
1615 // get base function gradient on rows
1616 auto t_row_base = row_data.getFTensor0N();
1617
1618 auto get_t_mat = [&](const int rr) {
1619 std::array<double *, FIELD_DIM * FIELD_DIM> ptrs;
1620 int s = 0;
1621 for (int i = 0; i != FIELD_DIM; ++i)
1622 for (int j = 0; j != FIELD_DIM; ++j, ++s)
1623 ptrs[s] = &(OpBase::locMat(rr + i, j));
1625 FIELD_DIM>(ptrs);
1626 };
1627
1628 auto t_u = getFTensor1FromMat<SPACE_DIM>(*uPtr);
1630
1634
1635 const double alpha_constant = alphaConstant();
1636 // loop over integration points
1637 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
1638 // take into account Jacobian
1639 const double alpha = t_w * vol * alpha_constant;
1640
1641 // loop over rows base functions
1642 int rr = 0;
1643 for (; rr != OpBase::nbRows / FIELD_DIM; ++rr) {
1644 // get matrix vec
1645 auto t_mat = get_t_mat(FIELD_DIM * rr);
1646 // get column base functions gradient at gauss point gg
1647 auto t_diff_col_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1648 // loop over columns
1649 for (int cc = 0; cc != OpBase::nbCols / FIELD_DIM; ++cc) {
1650 t_mat(I, L) +=
1651 alpha * t_row_base * t_kd(I, L) * (t_diff_col_base(k) * t_u(k));
1652 ++t_mat;
1653 ++t_diff_col_base;
1654 }
1655 ++t_row_base;
1656 }
1657 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
1658 ++t_row_base;
1659
1660 ++t_u;
1661 ++t_w; // move to another integration weight
1662 }
1664};
1665
1666} // namespace MoFEM
1667
1668#endif //__BILINEAR_FORMS_INTEGRATORS_HPP__
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 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 auto t_kd
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
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
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=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)
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=[](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=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)