386 {
388
389 auto get_sense_index = [
this]() {
return (
faceSense == 1) ? 0 : 1; };
390
394
395 int nb_integration_pts = OP::getGaussPts().size2();
396
398 auto t_kappa = getFTensor0FromVec<0>(*
kappaPtr);
400 auto t_face_normal = getFTensor1NormalsAtGaussPts();
401
402 int nb_base_functions = data.
getN().size2() / 3;
404 auto t_w = OP::getFTensor0IntegrationWeight();
405
406 auto next = [&]() {
407 ++t_P;
408 ++t_face_normal;
409 ++t_kappa;
410 ++t_delta_kappa;
411 ++t_w;
412 };
413
414 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
418 t_normalized_normal(
J) = t_normal(
J);
421 t_traction(
i) = t_P(
i,
J) * t_normalized_normal(
J);
422
423 auto fracture = [this](auto &t_traction, auto &t_normalized_normal,
424 auto &t_kappa, auto &t_delta_kappa, auto gc) {
425 double kappa =
static_cast<double>(t_kappa + t_delta_kappa);
432 t_gap_double(
i) = -t_gap(
i);
433 return t_gap_double;
434 };
435
437 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
438 int bb = 0;
439 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
440 double row_base = t_w * (t_row_base_fun(
J) * t_normal(
J));
441 t_nf(
i) += row_base * t_gap(
i);
442 ++t_nf;
443 ++t_row_base_fun;
444 }
445 for (; bb != nb_base_functions; ++bb)
446 ++t_row_base_fun;
447 };
448
449 assemble(fracture(t_traction, t_normalized_normal, t_kappa, t_delta_kappa,
451
452 next();
453 }
454
456 };
#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)
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.