476 {
478
479 auto get_sense_index = [
this]() {
return (
faceSense == 1) ? 0 : 1; };
480
485
486 int nb_integration_pts = OP::getGaussPts().size2();
487
489 auto t_kappa = getFTensor0FromVec<0>(*
kappaPtr);
491 auto t_face_normal = getFTensor1NormalsAtGaussPts();
492
493 int nb_base_functions = data.
getN().size2() / 3;
495 auto t_w = OP::getFTensor0IntegrationWeight();
496
497 auto next = [&]() {
498 ++t_P;
499 ++t_kappa;
500 ++t_delta_kappa;
501 ++t_face_normal;
502 ++t_w;
503 };
504
505 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
509 t_normalized_normal(
J) = t_normal(
J);
512 t_traction(
i) = t_P(
i,
J) * t_normalized_normal(
J);
513
514 auto fracture = [this](auto &t_traction, auto &t_normalized_normal,
515 auto &t_kappa, auto &t_delta_kappa, auto gc) {
516 double kappa =
static_cast<double>(t_kappa + t_delta_kappa);
524 t_dgap_double(
i,
j) = -t_dgap(
i,
j);
525 return t_dgap_double;
526 };
527
529 int rr = 0;
530 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
531 auto t_mat = getFTensor2FromArray<3, 3, 3>(OP::locMat, 3 * rr);
532 double row_base = t_w * (t_row_base_fun(
J) * t_normal(
J));
534 for (
int cc = 0; cc != nb_dofs /
SPACE_DIM; ++cc) {
535 double col_base = t_col_base_fun(
J) * t_normalized_normal(
J);
536 t_mat(
i,
j) += (row_base * col_base) * t_dgap(
i,
j);
537 ++t_mat;
538 ++t_col_base_fun;
539 }
540 ++t_row_base_fun;
541 }
542 for (; rr != nb_base_functions; ++rr)
543 ++t_row_base_fun;
544 };
545
546 assemble(fracture(t_traction, t_normalized_normal, t_kappa, t_delta_kappa,
548
549 next();
550 }
551
553 };
#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)
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.