435 {
437
438 auto get_sense_index = [
this]() {
return (
faceSense == 1) ? 0 : 1; };
439
443
444 int nb_integration_pts = OP::getGaussPts().size2();
445
447 auto t_kappa = getFTensor0FromVec<0>(*
kappaPtr);
449 auto t_face_normal = getFTensor1NormalsAtGaussPts();
450
451 int nb_base_functions = data.
getN().size2() / 3;
453 auto t_w = OP::getFTensor0IntegrationWeight();
454
455 auto next = [&]() {
456 ++t_P;
457 ++t_face_normal;
458 ++t_kappa;
459 ++t_delta_kappa;
460 ++t_w;
461 };
462
463 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
467 t_normalized_normal(
J) = t_normal(
J);
470 t_traction(
i) = t_P(
i,
J) * t_normalized_normal(
J);
471
472 auto fracture = [this](auto &t_traction, auto &t_normalized_normal,
473 auto &t_kappa, auto &t_delta_kappa, auto gc) {
474 double kappa =
static_cast<double>(t_kappa + t_delta_kappa);
479 true);
482 t_gap_double(
i) = -t_gap(
i) / 2.;
483 return t_gap_double;
484 };
485
487 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
488 int bb = 0;
489 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
490 double row_base = t_w * (t_row_base_fun(
J) * t_normal(
J));
491 t_nf(
i) += row_base * t_gap(
i);
492 ++t_nf;
493 ++t_row_base_fun;
494 }
495 for (; bb != nb_base_functions; ++bb)
496 ++t_row_base_fun;
497 };
498
499 assemble(fracture(t_traction, t_normalized_normal, t_kappa, t_delta_kappa,
501
502 next();
503 }
504
506 };
#define FTENSOR_INDEX(DIM, I)
Tensor1< T, Tensor_Dim > normalize()
#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()
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'J', DIM1 > J
static T getAlpha(const T &kappa, double gc, double min_stiffness)
static auto calculateGap(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
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....
auto getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
const VectorInt & getIndices() const
Get global indices of degrees of freedom on entity.