861 {
862
864
865 using CachePhi = boost::tuple<int, int, MatrixDouble>;
866
867 CGGUserPolynomialBase(boost::shared_ptr<CachePhi> cache_phi = nullptr);
868 ~CGGUserPolynomialBase() = default;
869
870 MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
871 BaseFunctionUnknownInterface **iface) const;
873 boost::shared_ptr<BaseFunctionCtx> ctx_ptr);
874
875private:
876 MatrixDouble shapeFun;
877 boost::shared_ptr<CachePhi> cachePhiPtr;
878
880};
881
882struct ContactTree;
883
891};
893
894using MatrixPtr = boost::shared_ptr<MatrixDouble>;
895using VectorPtr = boost::shared_ptr<VectorDouble>;
896
897using EntData = EntitiesFieldData::EntData;
900using FaceUserDataOperator = FaceElementForcesAndSourcesCore::UserDataOperator;
901using EleOnSide = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
902using SideEleOp = EleOnSide::UserDataOperator;
903
904struct PhysicalEquations;
905
906struct AnalyticalExprPython;
907
909 : public boost::enable_shared_from_this<DataAtIntegrationPts> {
910
911 MatrixDouble approxPAtPts;
912 MatrixDouble approxSigmaAtPts;
913 MatrixDouble divPAtPts;
914 MatrixDouble divSigmaAtPts;
915 MatrixDouble wL2AtPts;
916 MatrixDouble wL2DotAtPts;
917 MatrixDouble wL2DotDotAtPts;
918 MatrixDouble logStretchTensorAtPts;
919 MatrixDouble stretchTensorAtPts;
920 VectorDouble jacobianAtPts;
921 MatrixDouble tractionAtPts;
922
923 MatrixDouble diffStretchTensorAtPts;
924 VectorDouble detStretchTensorAtPts;
925 MatrixDouble detStretchTensorAtPts_du;
926 MatrixDouble logStretchDotTensorAtPts;
927 MatrixDouble gradLogStretchDotTensorAtPts;
928 MatrixDouble rotAxisAtPts;
929 MatrixDouble rotAxisDotAtPts;
930 MatrixDouble rotAxisGradDotAtPts;
931 MatrixDouble rotAxis0AtPts;
932 MatrixDouble WAtPts;
933 MatrixDouble W0AtPts;
934 MatrixDouble GAtPts;
935 MatrixDouble G0AtPts;
936 MatrixDouble wH1AtPts;
937 MatrixDouble XH1AtPts;
938 MatrixDouble contactL2AtPts;
939 MatrixDouble wGradH1AtPts;
940 MatrixDouble logStretch2H1AtPts;
941 MatrixDouble logStretchTotalTensorAtPts;
942
943 MatrixDouble hAtPts;
944 MatrixDouble hdOmegaAtPts;
945 MatrixDouble hdLogStretchAtPts;
946 MatrixDouble leviKirchhoffAtPts;
947 MatrixDouble leviKirchhoffdOmegaAtPts;
948 MatrixDouble leviKirchhoffdLogStreatchAtPts;
949 MatrixDouble leviKirchhoffPAtPts;
950 MatrixDouble adjointPdstretchAtPts;
951 MatrixDouble adjointPdUAtPts;
952 MatrixDouble adjointPdUdOmegaAtPts;
953 MatrixDouble adjointPdUdPAtPts;
954 MatrixDouble rotMatAtPts;
955 MatrixDouble PAtPts;
956 MatrixDouble SigmaAtPts;
957 VectorDouble energyAtPts;
958 MatrixDouble flowL2AtPts;
959
960 MatrixDouble varRotAxis;
961 MatrixDouble varLogStreach;
962 MatrixDouble varPiola;
963 MatrixDouble varDivPiola;
964 MatrixDouble varWL2;
965
966 MatrixDouble P_du;
967
968 MatrixDouble eigenVals;
969 MatrixDouble eigenVecs;
970 VectorInt nbUniq;
971 MatrixDouble eigenValsC;
972 MatrixDouble eigenVecsC;
973 VectorInt nbUniqC;
974
975 MatrixDouble matD;
976 MatrixDouble matAxiatorD;
977 MatrixDouble matDeviatorD;
978 MatrixDouble matInvD;
979
980 MatrixDouble internalStressAtPts;
981
982 MatrixDouble facePiolaAtPts;
983 MatrixDouble gradHybridDispAtPts;
984 MatrixDouble faceMaterialForceAtPts;
985 VectorDouble normalPressureAtPts;
986
989 double piolaScale = 1.;
990 double faceEnergy = 0.;
991
992 inline auto getPiolaScalePtr() {
993 return boost::shared_ptr<double>(shared_from_this(), &piolaScale);
994 }
995
997 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
998 &approxSigmaAtPts);
999 }
1001 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &approxPAtPts);
1002 }
1003
1005 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &divPAtPts);
1006 }
1007
1009 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &divSigmaAtPts);
1010 }
1011
1013 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &wL2AtPts);
1014 }
1015
1016 inline MatrixPtr getSmallWL2DotAtPts() {
1017 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &wL2DotAtPts);
1018 }
1019
1020 inline MatrixPtr getSmallWL2DotDotAtPts() {
1021 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &wL2DotDotAtPts);
1022 }
1023
1024 inline MatrixPtr getLogStretchTensorAtPts() {
1025 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
1026 &logStretchTensorAtPts);
1027 }
1028
1029 inline MatrixPtr getStretchTensorAtPts() {
1030 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
1031 &stretchTensorAtPts);
1032 }
1033
1034 inline MatrixPtr getLogStretchDotTensorAtPts() {
1035 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
1036 &logStretchDotTensorAtPts);
1037 }
1038
1039 inline MatrixPtr getGradLogStretchDotTensorAtPts() {
1040 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
1041 &gradLogStretchDotTensorAtPts);
1042 }
1043
1045 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &rotAxisAtPts);
1046 }
1047
1049 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &rotAxis0AtPts);
1050 }
1051
1053 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
1054 &rotAxisDotAtPts);
1055 }
1056
1057 inline MatrixPtr getRotAxisGradDotAtPts() {
1058 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
1059 &rotAxisGradDotAtPts);
1060 }
1061
1063 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &GAtPts);
1064 }
1065
1067 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &GAtPts);
1068 }
1069
1071 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matD);
1072 }
1073
1075 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matInvD);
1076 }
1077
1079 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matAxiatorD);
1080 }
1081
1083 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matDeviatorD);
1084 }
1085
1087 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &wH1AtPts);
1088 }
1089
1091 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &XH1AtPts);
1092 }
1093
1095 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &contactL2AtPts);
1096 }
1097
1098 inline MatrixPtr getSmallWGradH1AtPts() {
1099 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &wGradH1AtPts);
1100 }
1101
1103 return boost::shared_ptr<VectorDouble>(shared_from_this(), &jacobianAtPts);
1104 }
1105
1106 inline MatrixPtr getLeviKirchhoffAtPts() {
1107 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
1108 &leviKirchhoffAtPts);
1109 };
1110
1112 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varRotAxis);
1113 }
1114
1115 inline MatrixPtr getVarLogStreachPts() {
1116 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varLogStreach);
1117 }
1118
1120 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varPiola);
1121 }
1122
1124 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varDivPiola);
1125 }
1126
1128 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varWL2);
1129 }
1130
1131 inline MatrixPtr getInternalStressAtPts() {
1132 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
1133 &internalStressAtPts);
1134 }
1135
1137 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &facePiolaAtPts);
1138 }
1139
1140 inline MatrixPtr getGradHybridDispAtPts() {
1141 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
1142 &gradHybridDispAtPts);
1143 }
1144
1146 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &eigenVals);
1147 }
1148
1149 boost::shared_ptr<PhysicalEquations> physicsPtr;
1150};
1151
1153
1154
1155struct ExternalStrain;
1157struct PhysicalEquations {
1158
1165
1167
1170
1171 PhysicalEquations() = delete;
1172 PhysicalEquations(const int size_active, const int size_dependent)
1173 : activeVariables(size_active, 0),
1174 dependentVariablesPiola(size_dependent, 0),
1175 dependentVariablesPiolaDirevatives(size_dependent * size_active, 0) {}
1176 virtual ~PhysicalEquations() = default;
1177
1178 virtual MoFEMErrorCode recordTape(
const int tag, DTensor2Ptr *t_h) = 0;
1179
1181 returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs,
1182 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
1183 boost::shared_ptr<PhysicalEquations> physics_ptr);
1184
1186 returnOpSpatialPhysical(
const std::string &
field_name,
1187 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
1188 const double alpha_u);
1189
1191 returnOpSpatialPhysicalExternalStrain(
1193 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
1194 boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
1195 std::map<std::string, boost::shared_ptr<ScalingMethod>> smv);
1196
1198 std::string row_field, std::string col_field,
1199 boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha);
1200
1202 returnOpCalculateEnergy(boost::shared_ptr<DataAtIntegrationPts> data_ptr,
1203 boost::shared_ptr<double> total_energy_ptr);
1204
1206 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
1207 boost::shared_ptr<PhysicalEquations> physics_ptr);
1208
1210 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
1211 boost::shared_ptr<PhysicalEquations> physics_ptr);
1212
1214 returnOpSetScale(boost::shared_ptr<double> scale_ptr,
1215 boost::shared_ptr<PhysicalEquations> physics_ptr);
1216
1217 std::vector<double> activeVariables;
1218 std::vector<double> dependentVariablesPiola;
1219 std::vector<double> dependentVariablesPiolaDirevatives;
1220
1221
1222
1223
1224
1225 template <int S>
1226 inline static DTensor2Ptr get_VecTensor2(std::vector<double> &
v) {
1227 return DTensor2Ptr(&
v[S + 0], &
v[S + 1], &
v[S + 2], &
v[S + 3], &
v[S + 4],
1228 &
v[S + 5], &
v[S + 6], &
v[S + 7], &
v[S + 8]);
1229 }
1230
1231 template <int S>
1232 inline static DTensor0Ptr get_VecTensor0(std::vector<double> &
v) {
1233 return DTensor0Ptr(&
v[S + 0]);
1234 }
1235
1236 template <int S0>
1237 inline static DTensor3Ptr get_vecTensor3(std::vector<double> &
v,
1238 const int nba) {
1239
1240 const int A00 = nba * 0 + S0;
1241 const int A01 = nba * 1 + S0;
1242 const int A02 = nba * 2 + S0;
1243 const int A10 = nba * 3 + S0;
1244 const int A11 = nba * 4 + S0;
1245 const int A12 = nba * 5 + S0;
1246 const int A20 = nba * 6 + S0;
1247 const int A21 = nba * 7 + S0;
1248 const int A22 = nba * 8 + S0;
1249
1250 return DTensor3Ptr(
1251
1252 &
v[A00 + 0], &
v[A00 + 1], &
v[A00 + 2], &
v[A01 + 0], &
v[A01 + 1],
1253 &
v[A01 + 2], &
v[A02 + 0], &
v[A02 + 1], &
v[A02 + 2],
1254
1255 &
v[A10 + 0], &
v[A10 + 1], &
v[A10 + 2], &
v[A11 + 0], &
v[A11 + 1],
1256 &
v[A11 + 2], &
v[A12 + 0], &
v[A12 + 1], &
v[A12 + 2],
1257
1258 &
v[A20 + 0], &
v[A20 + 1], &
v[A20 + 2], &
v[A21 + 0], &
v[A21 + 1],
1259 &
v[A21 + 2], &
v[A22 + 0], &
v[A22 + 1], &
v[A22 + 2]
1260
1261 );
1262 }
1263
1264
1265
1266
1267
1268
1269
1270 inline DTensor2Ptr get_h() { return get_VecTensor2<0>(activeVariables); }
1271
1272
1273
1274
1275
1276
1277
1278 inline DTensor2Ptr get_P() {
1279 return get_VecTensor2<0>(dependentVariablesPiola);
1280 }
1281
1282
1283
1284
1285
1286
1287
1288 inline DTensor3Ptr get_P_dh0() {
1289 return get_vecTensor3<0>(dependentVariablesPiolaDirevatives,
1290 activeVariables.size());
1291 }
1292 inline DTensor3Ptr get_P_dh1() {
1293 return get_vecTensor3<3>(dependentVariablesPiolaDirevatives,
1294 activeVariables.size());
1295 }
1296 inline DTensor3Ptr get_P_dh2() {
1297 return get_vecTensor3<6>(dependentVariablesPiolaDirevatives,
1298 activeVariables.size());
1299 }
1300
1301
1302};
1303
1304struct BcDisp {
1305 BcDisp(std::string name, std::vector<double> attr,
Range faces);
1306 std::string blockName;
1308 VectorDouble3 vals;
1309 VectorInt3 flags;
1310};
1312
1313struct BcRot {
1314 BcRot(std::string name, std::vector<double> attr,
Range faces);
1315 std::string blockName;
1317 VectorDouble vals;
1318 double theta;
1319};
1320using BcRotVec = std::vector<BcRot>;
1321
1323
1324struct TractionBc {
1325 TractionBc(std::string name, std::vector<double> attr,
Range faces);
1326 std::string blockName;
1328 VectorDouble3 vals;
1329 VectorInt3 flags;
1330};
1332
1333struct NormalDisplacementBc {
1334 NormalDisplacementBc(std::string name, std::vector<double> attr,
Range faces);
1335 std::string blockName;
1337 double val;
1338};
1340
1341struct AnalyticalDisplacementBc {
1342 AnalyticalDisplacementBc(std::string name, std::vector<double> attr,
1344 std::string blockName;
1346 VectorInt3 flags;
1347};
1349
1350struct AnalyticalTractionBc {
1351 AnalyticalTractionBc(std::string name, std::vector<double> attr,
Range faces);
1352 std::string blockName;
1354 VectorInt3 flags;
1355};
1357
1358struct PressureBc {
1359 PressureBc(std::string name, std::vector<double> attr,
Range faces);
1360 std::string blockName;
1362 double val;
1363};
1365
1366struct ExternalStrain {
1367 ExternalStrain(std::string name, std::vector<double> attr,
Range ents);
1368 std::string blockName;
1370 double val;
1371 double bulkModulusK;
1372};
1373
1374 #ifdef ENABLE_PYTHON_BINDING
1375struct AnalyticalExprPython {
1376 AnalyticalExprPython() = default;
1377 virtual ~AnalyticalExprPython() = default;
1378
1380 MoFEMErrorCode evalAnalyticalDisp(
double delta_t,
double t, np::ndarray x,
1381 np::ndarray y, np::ndarray z,
1382 np::ndarray nx, np::ndarray ny,
1383 np::ndarray nz,
1384 const std::string &block_name,
1385 np::ndarray &analytical_expr);
1386 MoFEMErrorCode evalAnalyticalTraction(
double delta_t,
double t, np::ndarray x,
1387 np::ndarray y, np::ndarray z,
1388 np::ndarray nx, np::ndarray ny,
1389 np::ndarray nz,
1390 const std::string &block_name,
1391 np::ndarray &analytical_expr);
1392
1393 template <typename T>
1394 inline std::vector<T>
1395 py_list_to_std_vector(const boost::python::object &iterable) {
1396 return std::vector<T>(boost::python::stl_input_iterator<T>(iterable),
1397 boost::python::stl_input_iterator<T>());
1398 }
1399
1400private:
1401 bp::object mainNamespace;
1402 bp::object analyticalDispFun;
1403 bp::object analyticalTractionFun;
1404};
1405
1406extern boost::weak_ptr<AnalyticalExprPython> AnalyticalExprPythonWeakPtr;
1407 #endif
1408
1411
1412}
1413
1414#endif
1415
ForcesAndSourcesCore::UserDataOperator UserDataOperator
const double v
phase velocity of light in medium (cm/ns)
boost::shared_ptr< VectorDouble > VectorPtr
std::vector< AnalyticalTractionBc > AnalyticalTractionBcVec
std::vector< BcDisp > BcDispVec
std::vector< NormalDisplacementBc > NormalDisplacementBcVec
std::vector< BcRot > BcRotVec
std::vector< ExternalStrain > ExternalStrainVec
std::vector< AnalyticalDisplacementBc > AnalyticalDisplacementBcVec
std::vector< Range > TractionFreeBc
std::vector< PressureBc > PressureBcVec
std::vector< TractionBc > TractionBcVec
boost::shared_ptr< MatrixDouble > MatrixPtr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
constexpr double t
plate stiffness
constexpr auto field_name
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::FaceSideEle EleOnSide
Data on single entity (This is passed as argument to DataOperator::doWork)