995 {
997
1001
1002 const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
1003 auto &locMat = AssemblyBoundaryEleOp::locMat;
1004
1005 auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
1006 auto t_traction = getFTensor1FromMat<DIM>(
commonDataPtr->contactTraction);
1007 auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
1008
1009 auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
1010 auto t_row_base = row_data.getFTensor1N<3>();
1011 size_t nb_face_functions = row_data.getN().size2() / 3;
1012
1014 BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
1015 getFTensor1FromMat<DIM>(
commonDataPtr->contactDisp), nb_gauss_pts);
1017 BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
1018
1019 auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
1020
1021 auto ts_time = AssemblyBoundaryEleOp::getTStime();
1022 auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
1023
1024
1025 int block_id = 0;
1026
1027 auto v_sdf =
1029 m_spatial_coords, m_normals_at_pts, block_id);
1030
1031 auto m_grad_sdf =
1033 m_spatial_coords, m_normals_at_pts, block_id);
1034
1035 auto m_hess_sdf =
1037 m_spatial_coords, m_normals_at_pts, block_id);
1038
1039 auto t_sdf = getFTensor0FromVec(v_sdf);
1040 auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
1041 auto t_hess_sdf = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
1042
1043 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1044
1045 double jacobian = 1.;
1047 jacobian = 2. * M_PI * t_coords(0);
1048 }
1049 const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
1050
1051 auto tn = -t_traction(
i) * t_grad_sdf(
i);
1053
1055 t_cP(
i,
j) = (
c * t_grad_sdf(
i)) * t_grad_sdf(
j);
1058
1061
1065 (t_hess_sdf(
i,
j) * (t_grad_sdf(
k) * t_traction(
k)) +
1066 t_grad_sdf(
i) * t_hess_sdf(
k,
j) * t_traction(
k)) +
1067 c * t_sdf * t_hess_sdf(
i,
j);
1068 }
1069
1070 size_t rr = 0;
1071 for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
1072
1073 auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * rr);
1074
1075 const double row_base = t_row_base(
i) * t_normal(
i);
1076
1077 auto t_col_base = col_data.getFTensor0N(gg, 0);
1078 for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / DIM; ++cc) {
1079 const double beta = alpha * row_base * t_col_base;
1080
1081 t_mat(
i,
j) -= beta * t_res_dU(
i,
j);
1082
1083 ++t_col_base;
1084 ++t_mat;
1085 }
1086
1087 ++t_row_base;
1088 }
1089 for (; rr < nb_face_functions; ++rr)
1090 ++t_row_base;
1091
1092 ++t_traction;
1093 ++t_coords;
1094 ++t_w;
1095 ++t_normal;
1096 ++t_sdf;
1097 ++t_grad_sdf;
1098 ++t_hess_sdf;
1099 }
1100
1102}
#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
const double c
speed of light (cm/ns)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.