511 {
513
514
516 const auto in_the_loop =
518
519
522
523
524 auto t_normal = getFTensor1Normal();
525 t_normal.normalize();
526
527
528
530
531
532 const size_t nb_integration_pts = getGaussPts().size2();
533
534
535
536 const double beta =
static_cast<double>(
nitsche) / (in_the_loop + 1);
537
538 auto integrate = [&](auto sense_row, auto &row_ind, auto &row_diff,
539 auto &row_diff2, auto sense_col, auto &col_ind,
540 auto &col_diff, auto &col_diff2) {
542
543
544
545 const auto nb_rows = row_ind.size();
546 const auto nb_cols = col_ind.size();
547
548 const auto nb_row_base_functions = row_diff.size2() /
SPACE_DIM;
549
550 if (nb_cols) {
551
552
553 locMat.resize(nb_rows, nb_cols,
false);
555
556
559
560 auto t_w = getFTensor0IntegrationWeight();
561
562
563 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
564
565
566
567 const double alpha = getMeasure() * t_w;
568 auto t_mat =
locMat.data().begin();
569
570
571 size_t rr = 0;
572 for (; rr != nb_rows; ++rr) {
573
575 t_mv(
i,
j) = t_D(
i,
j,
k,
l) * t_diff2_row_base(
k,
l);
576
577
579 t_vn_plus(
i,
j) = beta * (
phi * t_mv(
i,
j) / p);
581 t_vn(
i,
j) = (t_diff_row_base(
j) * (t_normal(
i) * sense_row)) -
583
584
587
588
589 for (size_t cc = 0; cc != nb_cols; ++cc) {
590
592 t_mu(
i,
j) = t_D(
i,
j,
k,
l) * t_diff2_col_base(
k,
l);
593
594
596 t_un(
i,
j) = -p * ((t_diff_col_base(
j) * (t_normal(
i) * sense_col) -
597 beta * t_mu(
i,
j) / p));
598
599
600 *t_mat -= alpha * (t_vn(
i,
j) * t_un(
i,
j));
601 *t_mat -= alpha * (t_vn_plus(
i,
j) * (beta * t_mu(
i,
j)));
602
603
604 ++t_diff_col_base;
605 ++t_diff2_col_base;
606 ++t_mat;
607 }
608
609
610 ++t_diff_row_base;
611 ++t_diff2_row_base;
612 }
613
614
615
616
617
618
619 for (; rr < nb_row_base_functions; ++rr) {
620 ++t_diff_row_base;
621 ++t_diff2_row_base;
622 }
623
624 ++t_w;
625 }
626
627
628 CHKERR ::MatSetValues(getKSPB(), nb_rows, &*row_ind.begin(),
629 col_ind.size(), &*col_ind.begin(),
630 &*
locMat.data().begin(), ADD_VALUES);
631 }
632
634 };
635
636
638
639 const auto sense_row =
senseMap[s0];
640
642
644 const auto sense_col =
senseMap[s1];
645
647
650
653
654 );
655 }
656 }
657 }
658 }
659
661}
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
static FTensor::Ddg< FTensor::PackPtr< T *, 1 >, Tensor_Dim01, Tensor_Dim23 > getFTensor4DdgFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 4 on first two and last indices from form data matrix.
auto get_diff_ntensor(T &base_mat)
std::array< std::vector< MatrixDouble >, 2 > diff2BaseSideMap
FTensor::Index< 'j', SPACE_DIM > j
std::array< std::vector< MatrixDouble >, 2 > diffBaseSideMap
FTensor::Index< 'k', SPACE_DIM > k
FTensor::Index< 'i', SPACE_DIM > i
constexpr int SPACE_DIM
dimension of space
FTensor::Index< 'l', SPACE_DIM > l
std::array< double, 2 > areaMap
std::array< std::vector< VectorInt >, 2 > indicesSideMap
indices on rows for left hand-side
std::array< int, 2 > senseMap
auto get_diff2_ntensor(T &base_mat)
MatrixDouble locMat
local operator matrix