494 {
496
497 auto get_sense_index = [
this]() {
return (
faceSense == 1) ? 0 : 1; };
498
503
504 int nb_integration_pts = OP::getGaussPts().size2();
505
507 auto t_kappa = getFTensor0FromVec<0>(*
kappaPtr);
509 auto t_face_normal = getFTensor1NormalsAtGaussPts();
510
511 int nb_base_functions = data.
getN().size2() / 3;
513 auto t_w = OP::getFTensor0IntegrationWeight();
514
515 auto next = [&]() {
516 ++t_P;
517 ++t_kappa;
518 ++t_delta_kappa;
519 ++t_face_normal;
520 ++t_w;
521 };
522
523 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
527 t_normalized_normal(
J) = t_normal(
J);
530 t_traction(
i) = t_P(
i,
J) * t_normalized_normal(
J);
531
532 auto fracture = [this](auto &t_traction, auto &t_normalized_normal,
533 auto &t_kappa, auto &t_delta_kappa, auto gc) {
534 double kappa =
static_cast<double>(t_kappa + t_delta_kappa);
539 true);
543 t_dgap_double(
i,
j) = -t_dgap(
i,
j);
544 return t_dgap_double;
545 };
546
548 int rr = 0;
549 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
550 auto t_mat = getFTensor2FromArray<3, 3, 3>(OP::locMat, 3 * rr);
551 double row_base = t_w * (t_row_base_fun(
J) * t_normal(
J));
553 for (
int cc = 0; cc != nb_dofs /
SPACE_DIM; ++cc) {
554 double col_base = t_col_base_fun(
J) * t_normalized_normal(
J);
555 t_mat(
i,
j) += (row_base * col_base) * t_dgap(
i,
j);
556 ++t_mat;
557 ++t_col_base_fun;
558 }
559 ++t_row_base_fun;
560 }
561 for (; rr != nb_base_functions; ++rr)
562 ++t_row_base_fun;
563 };
564
565 assemble(fracture(t_traction, t_normalized_normal, t_kappa, t_delta_kappa,
567
568 next();
569 }
570
572 };
#define FTENSOR_INDEX(DIM, I)
Tensor1< T, Tensor_Dim > normalize()
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'J', DIM1 > J
FTensor::Index< 'j', 3 > j
static T getAlpha(const T &kappa, double gc, double min_stiffness)
static auto calculateDiffGapDTraction(const FTensor::Tensor1< T, 3 > &t_traction, FTensor::Tensor1< double, 3 > &t_n_normalize, double alpha, double beta, bool sign_sensitive=true)
boost::shared_ptr< MatrixDouble > fluxMatPtr
boost::shared_ptr< double > gcPtr
boost::shared_ptr< VectorDouble > kappaDeltaPtr
boost::shared_ptr< VectorDouble > kappaPtr
static double min_stiffness
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 degrees of freedom on entity.