508 {
510
511
513 const auto in_the_loop =
515
516
519
520
521 auto t_normal = getFTensor1Normal();
522 t_normal.normalize();
523
524
525
526 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*
dMatPtr);
527
528
529 const size_t nb_integration_pts = getGaussPts().size2();
530
531
532
533 const double beta =
static_cast<double>(
nitsche) / (in_the_loop + 1);
534
535 auto integrate = [&](auto sense_row, auto &row_ind, auto &row_diff,
536 auto &row_diff2, auto sense_col, auto &col_ind,
537 auto &col_diff, auto &col_diff2) {
539
540
541
542 const auto nb_rows = row_ind.size();
543 const auto nb_cols = col_ind.size();
544
545 const auto nb_row_base_functions = row_diff.size2() /
SPACE_DIM;
546
547 if (nb_cols) {
548
549
550 locMat.resize(nb_rows, nb_cols,
false);
552
553
556
557 auto t_w = getFTensor0IntegrationWeight();
558
559
560 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
561
562
563
564 const double alpha = getMeasure() * t_w;
565 auto t_mat =
locMat.data().begin();
566
567
568 size_t rr = 0;
569 for (; rr != nb_rows; ++rr) {
570
572 t_mv(
i,
j) = t_D(
i,
j,
k,
l) * t_diff2_row_base(
k,
l);
573
574
576 t_vn_plus(
i,
j) = beta * (
phi * t_mv(
i,
j) / p);
578 t_vn(
i,
j) = (t_diff_row_base(
j) * (t_normal(
i) * sense_row)) -
580
581
584
585
586 for (size_t cc = 0; cc != nb_cols; ++cc) {
587
589 t_mu(
i,
j) = t_D(
i,
j,
k,
l) * t_diff2_col_base(
k,
l);
590
591
593 t_un(
i,
j) = -p * ((t_diff_col_base(
j) * (t_normal(
i) * sense_col) -
594 beta * t_mu(
i,
j) / p));
595
596
597 *t_mat -= alpha * (t_vn(
i,
j) * t_un(
i,
j));
598 *t_mat -= alpha * (t_vn_plus(
i,
j) * (beta * t_mu(
i,
j)));
599
600
601 ++t_diff_col_base;
602 ++t_diff2_col_base;
603 ++t_mat;
604 }
605
606
607 ++t_diff_row_base;
608 ++t_diff2_row_base;
609 }
610
611
612
613
614
615
616 for (; rr < nb_row_base_functions; ++rr) {
617 ++t_diff_row_base;
618 ++t_diff2_row_base;
619 }
620
621 ++t_w;
622 }
623
624
625 CHKERR ::MatSetValues(getKSPB(), nb_rows, &*row_ind.begin(),
626 col_ind.size(), &*col_ind.begin(),
627 &*
locMat.data().begin(), ADD_VALUES);
628 }
629
631 };
632
633
635
636 const auto sense_row =
senseMap[s0];
637
639
641 const auto sense_col =
senseMap[s1];
642
644
647
650
651 );
652 }
653 }
654 }
655 }
656
658}
#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 ...
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