403 {
405
406 auto get_sense_index = [
this]() {
return (
faceSense == 1) ? 0 : 1; };
407
411
412 int nb_integration_pts = OP::getGaussPts().size2();
413
415 auto t_kappa = getFTensor0FromVec<0>(*
kappaPtr);
417 auto t_face_normal = getFTensor1NormalsAtGaussPts();
418
419 int nb_base_functions = data.
getN().size2() / 3;
421 auto t_w = OP::getFTensor0IntegrationWeight();
422
423 auto next = [&]() {
424 ++t_P;
425 ++t_face_normal;
426 ++t_kappa;
427 ++t_delta_kappa;
428 ++t_w;
429 };
430
431 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
435 t_normalized_normal(
J) = t_normal(
J);
438 t_traction(
i) = t_P(
i,
J) * t_normalized_normal(
J);
439
440 auto fracture = [this](auto &t_traction, auto &t_normalized_normal,
441 auto &t_kappa, auto &t_delta_kappa, auto gc) {
442 double kappa =
static_cast<double>(t_kappa + t_delta_kappa);
447 true);
450 t_gap_double(
i) = -t_gap(
i);
451 return t_gap_double;
452 };
453
455 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
456 int bb = 0;
457 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
458 double row_base = t_w * (t_row_base_fun(
J) * t_normal(
J));
459 t_nf(
i) += row_base * t_gap(
i);
460 ++t_nf;
461 ++t_row_base_fun;
462 }
463 for (; bb != nb_base_functions; ++bb)
464 ++t_row_base_fun;
465 };
466
467 assemble(fracture(t_traction, t_normalized_normal, t_kappa, t_delta_kappa,
469
470 next();
471 }
472
474 };
#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....
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.