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\).
3336 {
3338
3351
3352 dataAtPts->faceMaterialForceAtPts.resize(3, getGaussPts().size2(),
false);
3353 dataAtPts->normalPressureAtPts.resize(getGaussPts().size2(),
false);
3354 if (getNinTheLoop() == 0) {
3355 dataAtPts->faceMaterialForceAtPts.clear();
3357 }
3358 auto loop_size = getLoopSize();
3359 if (loop_size == 1) {
3360 auto numebered_fe_ptr = getSidePtrFE()->numeredEntFiniteElementPtr;
3361 auto pstatus = numebered_fe_ptr->getPStatus();
3362 if (pstatus & (PSTATUS_SHARED | PSTATUS_MULTISHARED)) {
3363 loop_size = 2;
3364 }
3365 }
3366
3368
3369 auto t_normal = getFTensor1NormalsAtGaussPts();
3370 auto t_T = getFTensor1FromMat<SPACE_DIM>(
3372 auto t_p =
3374 auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
dataAtPts->approxPAtPts);
3375
3376
3377 auto t_u_gamma =
3378 getFTensor1FromMat<SPACE_DIM>(
dataAtPts->hybridDispAtPts);
3379 auto t_grad_u_gamma =
3380 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
dataAtPts->gradHybridDispAtPts);
3381 auto t_strain =
3382 getFTensor2SymmetricFromMat<SPACE_DIM>(
dataAtPts->logStretchTensorAtPts);
3383 auto t_omega = getFTensor1FromMat<3>(
dataAtPts->rotAxisAtPts);
3384
3390
3391 auto next = [&]() {
3392 ++t_normal;
3393 ++t_P;
3394
3395 ++t_omega;
3396 ++t_u_gamma;
3397 ++t_grad_u_gamma;
3398 ++t_strain;
3399 ++t_T;
3400 ++t_p;
3401 };
3402
3405 for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3406 t_N(
I) = t_normal(
I);
3408
3411 t_grad_u(
i,
j) = t_R(
i,
j) + t_strain(
i,
j);
3412
3414 t_N(
J) * (t_grad_u(
i,
I) * t_P(
i,
J)) / loop_size;
3415
3416
3417 t_T(
I) -= t_N(
I) * ((t_strain(
i,
K) * t_P(
i,
K)) / 2.) / loop_size;
3418
3419 t_p += t_N(
I) * (t_N(
J) * (t_grad_u_gamma(
i,
I) * t_P(
i,
J))) / loop_size;
3420
3421 next();
3422 }
3423 break;
3425 for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3426
3427
3428 t_N(
I) = t_normal(
I);
3430
3431
3433 t_strain(
i,
j) - 0.5 * (t_grad_u_gamma(
i,
j) + t_grad_u_gamma(
j,
i));
3434
3435
3436 t_grad_u(
i,
J) = t_grad_u_gamma(
i,
J) + (2 * t_R(
i,
K) * t_N(
K) -
3437 (t_R(
k,
L) * t_N(
k) * t_N(
L)) * t_N(
i)) *
3439
3440 t_T(
I) += t_N(
J) * (t_grad_u(
i,
I) * t_P(
i,
J)) / loop_size;
3441
3442
3443 t_T(
I) -= t_N(
I) * ((t_strain(
i,
K) * t_P(
i,
K)) / 2.) / loop_size;
3444
3445
3446 t_p += t_N(
I) * (t_N(
J) * (t_grad_u_gamma(
i,
I) * t_P(
i,
J))) / loop_size;
3447
3448 next();
3449 }
3450 break;
3451
3452 default:
3454 "Grffith energy release "
3455 "selector not implemented");
3456 };
3457
3458#ifndef NDEBUG
3459 auto side_fe_ptr = getSidePtrFE();
3460 auto side_fe_mi_ptr = side_fe_ptr->numeredEntFiniteElementPtr;
3461 auto pstatus = side_fe_mi_ptr->getPStatus();
3462 if (pstatus) {
3463 auto owner = side_fe_mi_ptr->getOwnerProc();
3465 << "OpFaceSideMaterialForce: owner proc is not 0, owner proc: " << owner
3466 << " " << getPtrFE()->mField.get_comm_rank() << " n in the loop "
3467 << getNinTheLoop() << " loop size " << getLoopSize();
3468 }
3469#endif
3470
3472}
#define FTENSOR_INDEX(DIM, I)
constexpr int SPACE_DIM
[Define dimension]
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