Caluclate face material force and normal pressure at gauss points.
Procedure (pointwise on S): 1) Normalize the normal \(\mathbf n=\mathbf N/\|\mathbf N\|\). 2) Form the residual \(R=\varepsilon-\operatorname{sym}(u^\Gamma)\), where \(\operatorname{sym}(A)=\tfrac12(A+A^\top)\). 3) Recover the normal directional derivative (a vector) \(\mathbf v=\partial_{\mathbf n}u=2R\mathbf n-(\mathbf n^\top R\,\mathbf n)\,\mathbf n\). 4) Assemble the full gradient \(U = u^\Gamma + \mathbf v\otimes \mathbf n\).
3368 {
3370
3383
3384 dataAtPts->faceMaterialForceAtPts.resize(3, getGaussPts().size2(),
false);
3385 dataAtPts->normalPressureAtPts.resize(getGaussPts().size2(),
false);
3386 if (getNinTheLoop() == 0) {
3387 dataAtPts->faceMaterialForceAtPts.clear();
3389 }
3390 auto loop_size = getLoopSize();
3391 if (loop_size == 1) {
3392 auto numebered_fe_ptr = getSidePtrFE()->numeredEntFiniteElementPtr;
3393 auto pstatus = numebered_fe_ptr->getPStatus();
3394 if (pstatus & (PSTATUS_SHARED | PSTATUS_MULTISHARED)) {
3395 loop_size = 2;
3396 }
3397 }
3398
3400
3401 auto t_normal = getFTensor1NormalsAtGaussPts();
3402 auto t_T = getFTensor1FromMat<SPACE_DIM>(
3404 auto t_p =
3406 auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
dataAtPts->approxPAtPts);
3407
3408
3409 auto t_u_gamma =
3410 getFTensor1FromMat<SPACE_DIM>(
dataAtPts->hybridDispAtPts);
3411 auto t_grad_u_gamma =
3412 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
dataAtPts->gradHybridDispAtPts);
3413 auto t_strain =
3414 getFTensor2SymmetricFromMat<SPACE_DIM>(
dataAtPts->logStretchTensorAtPts);
3415 auto t_omega = getFTensor1FromMat<3>(
dataAtPts->rotAxisAtPts);
3416
3422
3423 auto next = [&]() {
3424 ++t_normal;
3425 ++t_P;
3426
3427 ++t_omega;
3428 ++t_u_gamma;
3429 ++t_grad_u_gamma;
3430 ++t_strain;
3431 ++t_T;
3432 ++t_p;
3433 };
3434
3437 for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3438 t_N(
I) = t_normal(
I);
3440
3443 t_grad_u(
i,
j) = t_R(
i,
j) + t_strain(
i,
j);
3444
3446 t_N(
J) * (t_grad_u(
i,
I) * t_P(
i,
J)) / loop_size;
3447
3448
3449 t_T(
I) -= t_N(
I) * ((t_strain(
i,
K) * t_P(
i,
K)) / 2.) / loop_size;
3450
3451 t_p += t_N(
I) * (t_N(
J) * (t_grad_u_gamma(
i,
I) * t_P(
i,
J))) / loop_size;
3452
3453 next();
3454 }
3455 break;
3457 for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3458
3459
3460 t_N(
I) = t_normal(
I);
3462
3463
3465 t_strain(
i,
j) - 0.5 * (t_grad_u_gamma(
i,
j) + t_grad_u_gamma(
j,
i));
3466
3467
3469 t_grad_u_gamma(
i,
J) +
3470 (2 * t_R(
i,
K) * t_N(
K) - (t_R(
k,
L) * t_N(
k) * t_N(
L)) * t_N(
i)) *
3472
3473 t_T(
I) += t_N(
J) * (t_grad_u(
i,
I) * t_P(
i,
J)) / loop_size;
3474
3475
3476 t_T(
I) -= t_N(
I) * ((t_strain(
i,
K) * t_P(
i,
K)) / 2.) / loop_size;
3477
3478
3479 t_p += t_N(
I) * (t_N(
J) * (t_grad_u_gamma(
i,
I) * t_P(
i,
J))) / loop_size;
3480
3481 next();
3482 }
3483 break;
3484
3485 default:
3487 "Grffith energy release "
3488 "selector not implemented");
3489 };
3490
3491#ifndef NDEBUG
3492 auto side_fe_ptr = getSidePtrFE();
3493 auto side_fe_mi_ptr = side_fe_ptr->numeredEntFiniteElementPtr;
3494 auto pstatus = side_fe_mi_ptr->getPStatus();
3495 if (pstatus) {
3496 auto owner = side_fe_mi_ptr->getOwnerProc();
3498 << "OpFaceSideMaterialForce: owner proc is not 0, owner proc: " << owner
3499 << " " << getPtrFE()->mField.get_comm_rank() << " n in the loop "
3500 << getNinTheLoop() << " loop size " << getLoopSize();
3501 }
3502#endif
3503
3505}
#define FTENSOR_INDEX(DIM, I)
Kronecker Delta class symmetric.
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()
#define MOFEM_LOG(channel, severity)
Log.
FTensor::Index< 'i', SPACE_DIM > i
const double n
refractive index of diffusive medium
FTensor::Index< 'J', DIM1 > J
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
constexpr IntegrationType I
FTensor::Index< 'm', 3 > m
static enum EnergyReleaseSelector energyReleaseSelector