1034 {
1037
1038#ifndef NDEBUG
1041 "Topological data pointer is null");
1044 "Broken displacement data pointer is null");
1047 "Hybrid displacement pointer is null");
1050 "Adjoint hybrid displacement pointer is null");
1051#endif
1052
1053 const int nb_dofs =
data.getIndices().size();
1054 if (!nb_dofs)
1056
1057 const int nb_integration_pts = getGaussPts().size2();
1058 const int nb_base_functions =
data.getN().size2();
1059
1060#ifndef NDEBUG
1061 if (this->
nF.size() != nb_dofs)
1063 "Size of nF %ld != nb_dofs %d", this->
nF.size(), nb_dofs);
1064 if (
data.getDiffN().size1() != nb_integration_pts)
1066 "Differential of base functions should have the same number of "
1067 "integration points as the data");
1068 if (
data.getDiffN().size2() != nb_base_functions * 2)
1070 "Differential of base functions should have the same number of "
1071 "base functions as the data");
1072#endif
1073
1077
1078 auto &coords = getCoords();
1079
1080
1082 std::get<2>(Tools::getTricircumcenter3d(coords.data().data()));
1083
1084 auto t_w = getFTensor0IntegrationWeight();
1085 auto t_tangent1 = getFTensor1Tangent1AtGaussPts();
1086 auto t_tangent2 = getFTensor1Tangent2AtGaussPts();
1088 auto t_var_u_hybrid =
1090 auto t_diff_base =
data.getFTensor1DiffN<2>();
1091
1092 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1096
1097 double area = std::sqrt(t_normal(
i) * t_normal(
i));
1099 t_da(
i) = t_normal(
i) / area;
1100 area /= 2.;
1102
1103 double tau_density = 0;
1106 bd.getFlux(), nb_integration_pts);
1108 bd.getVarFlux(), nb_integration_pts);
1109 for (int ss = 0; ss != gg; ++ss) {
1110 ++t_u_broken;
1111 ++t_var_u_broken;
1112 }
1113
1114
1115
1116
1117 const double hybrid_hybrid = t_var_u_hybrid(
i) * t_u_hybrid(
i);
1118 const double broken_broken = t_var_u_broken(
i) * t_u_broken(
i);
1119 const double hybrid_broken = -t_var_u_hybrid(
i) * t_u_broken(
i);
1120 const double broken_hybrid = -t_var_u_broken(
i) * t_u_hybrid(
i);
1121 tau_density +=
1122 hybrid_hybrid + broken_broken + hybrid_broken + broken_hybrid;
1123 }
1124
1126 locJ += t_w * tau * area * tau_density;
1127
1128 auto t_nf = getFTensor1FromArray<SPACE_DIM, SPACE_DIM>(
nF);
1129 int rr = 0;
1130 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
1135
1137 t_w *
alphaTau * tau_density * (t_da(
i) * t_normal_dX(
i,
I)) /
h;
1138 ++t_diff_base;
1139 ++t_nf;
1140 }
1141 for (; rr != nb_base_functions; ++rr)
1142 ++t_diff_base;
1143
1144 ++t_w;
1145 ++t_tangent1;
1146 ++t_tangent2;
1147 ++t_u_hybrid;
1148 ++t_var_u_hybrid;
1149 }
1150
1152}
boost::shared_ptr< DataAtIntegrationPts > data
#define FTENSOR_INDEXES(DIM,...)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#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', 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
auto getFTensor1FromMat(M &data, int rr=0, int cc=0)
Get tensor rank 1 (vector) form data matrix.
constexpr IntegrationType I
boost::shared_ptr< TopologicalData > topoData
VectorDouble nF
local right hand side vector
Public Member Functions inherited from EshelbianPlasticity::OpAssembleTopologicalObjectiveDerivativeImplBase< OpAssembleFace >