Caluclate face material force and normal pressure at gauss points.
Reconstruct the full gradient \(U=\nabla u\) on a surface from the symmetric part and the surface gradient.
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\).
3359 {
3361
3374
3375 dataAtPts->faceMaterialForceAtPts.resize(3, getGaussPts().size2(),
false);
3376 dataAtPts->normalPressureAtPts.resize(getGaussPts().size2(),
false);
3377 if (getNinTheLoop() == 0) {
3378 dataAtPts->faceMaterialForceAtPts.clear();
3380 }
3381 auto loop_size = getLoopSize();
3382 if (loop_size == 1) {
3383 auto numebered_fe_ptr = getSidePtrFE()->numeredEntFiniteElementPtr;
3384 auto pstatus = numebered_fe_ptr->getPStatus();
3385 if (pstatus & (PSTATUS_SHARED | PSTATUS_MULTISHARED)) {
3386 loop_size = 2;
3387 }
3388 }
3389
3391
3392 auto t_normal = getFTensor1NormalsAtGaussPts();
3393 auto t_T = getFTensor1FromMat<SPACE_DIM>(
3395 auto t_p =
3397 auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
dataAtPts->approxPAtPts);
3398
3399
3400 auto t_u_gamma =
3401 getFTensor1FromMat<SPACE_DIM>(
dataAtPts->hybridDispAtPts);
3402 auto t_grad_u_gamma =
3403 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
dataAtPts->gradHybridDispAtPts);
3404 auto t_strain =
3405 getFTensor2SymmetricFromMat<SPACE_DIM>(
dataAtPts->logStretchTensorAtPts);
3406 auto t_omega = getFTensor1FromMat<3>(
dataAtPts->rotAxisAtPts);
3407
3413
3414 auto next = [&]() {
3415 ++t_normal;
3416 ++t_P;
3417
3418 ++t_omega;
3419 ++t_u_gamma;
3420 ++t_grad_u_gamma;
3421 ++t_strain;
3422 ++t_T;
3423 ++t_p;
3424 };
3425
3428 for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3429 t_N(
I) = t_normal(
I);
3431
3434 t_grad_u(
i,
j) = t_R(
i,
j) + t_strain(
i,
j);
3435
3437 t_N(
J) * (t_grad_u(
i,
I) * t_P(
i,
J)) / loop_size;
3438
3439
3440 t_T(
I) -= t_N(
I) * ((t_strain(
i,
K) * t_P(
i,
K)) / 2.) / loop_size;
3441
3443 (t_N(
J) * ((
t_kd(
i,
I) + t_grad_u_gamma(
i,
I)) * t_P(
i,
J))) /
3444 loop_size;
3445
3446 next();
3447 }
3448 break;
3450 for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3451
3452
3453 t_N(
I) = t_normal(
I);
3455
3456
3458 t_strain(
i,
j) - 0.5 * (t_grad_u_gamma(
i,
j) + t_grad_u_gamma(
j,
i));
3459
3460
3461 t_grad_u(
i,
J) = t_grad_u_gamma(
i,
J) + (2 * t_R(
i,
K) * t_N(
K) -
3462 (t_R(
k,
L) * t_N(
k) * t_N(
L)) * t_N(
i)) *
3464
3465 t_T(
I) += t_N(
J) * (t_grad_u(
i,
I) * t_P(
i,
J)) / loop_size;
3466
3467
3468 t_T(
I) -= t_N(
I) * ((t_strain(
i,
K) * t_P(
i,
K)) / 2.) / loop_size;
3469
3470
3472 (t_N(
J) * ((
t_kd(
i,
I) + t_grad_u_gamma(
i,
I)) * t_P(
i,
J))) /
3473 loop_size;
3474
3475 next();
3476 }
3477 break;
3478
3479 default:
3481 "Grffith energy release "
3482 "selector not implemented");
3483 };
3484
3485#ifndef NDEBUG
3486 auto side_fe_ptr = getSidePtrFE();
3487 auto side_fe_mi_ptr = side_fe_ptr->numeredEntFiniteElementPtr;
3488 auto pstatus = side_fe_mi_ptr->getPStatus();
3489 if (pstatus) {
3490 auto owner = side_fe_mi_ptr->getOwnerProc();
3492 << "OpFaceSideMaterialForce: owner proc is not 0, owner proc: " << owner
3493 << " " << getPtrFE()->mField.get_comm_rank() << " n in the loop "
3494 << getNinTheLoop() << " loop size " << getLoopSize();
3495 }
3496#endif
3497
3499}
#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