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