Operator for linear form, usually to calculate values on right hand side.
965 {
967 if (row_type != MBTET && row_type != MBHEX)
970
975
976 if (SEpEp) {
977
979 auto &col_name, auto row_side, auto col_side,
980 auto row_type, auto col_type) {
981 return add_bmc.get<0>().find(boost::make_tuple(
982 row_name, col_name, row_type, col_type, row_side, col_side));
983 };
984
986 auto &col_name, auto row_side, auto col_side,
987 auto row_type,
auto col_type,
const auto &
m,
988 const auto &row_ind, const auto &col_ind) {
989 auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
990 row_type, col_type);
991 if (it != add_bmc.get<0>().end()) {
993 it->setInd(row_ind, col_ind);
994 it->setSetAtElement();
995 return it;
996 } else {
997 auto p_it = add_bmc.insert(
BlockMatData(row_name, col_name, row_type,
998 col_type, row_side, col_side,
999 m, row_ind, col_ind));
1000 return p_it.first;
1001 }
1002 };
1003
1005 auto &col_name, auto row_side, auto col_side,
1006 auto row_type,
auto col_type,
const auto &
m,
1007 const auto &row_ind, const auto &col_ind) {
1008 auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
1009 row_type, col_type);
1010 if (it != add_bmc.get<0>().end()) {
1012 it->setInd(row_ind, col_ind);
1013 it->setSetAtElement();
1014 return it;
1015 } else {
1016 auto p_it = add_bmc.insert(
BlockMatData(row_name, col_name, row_type,
1017 col_type, row_side, col_side,
1018 m, row_ind, col_ind));
1019 return p_it.first;
1020 }
1021 };
1022
1023 auto assemble_block = [&](
auto &
bit, Mat S) {
1028
1030 &*cind.begin(), &*
m.data().begin(), ADD_VALUES);
1031
1033 };
1034
1037 const int nb =
m.size1();
1038
1039 inv.resize(nb, nb, false);
1040 inv.clear();
1041 for (int cc = 0; cc != nb; ++cc)
1042 inv(cc, cc) = -1;
1043
1044 iPIV.resize(nb,
false);
1046 const auto info =
lapack_dsysv(
'L', nb, nb, &*
m.data().begin(), nb,
1047 &*
iPIV.begin(), &*inv.data().begin(), nb,
1049
1050
1051
1052
1054 };
1055
1057 auto &inv) {
1059
1061 bmc.get<1>().find(boost::make_tuple(field, field, row_type, row_type));
1062 if (
bit != bmc.get<1>().end()) {
1063
1065 CHKERR invert_symm_mat(
m, inv);
1066
1067 } else
1069 "%s matrix not found", field.c_str());
1070
1072 };
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1127 std::string field, AO ao,
1130
1131 for (
auto &
bit : add_bmc) {
1132 bit.unSetAtElement();
1134 }
1135
1136 for (
auto &
bit : bmc) {
1137 if (
bit.setAtElement &&
bit.rowField != field &&
1138 bit.colField != field) {
1142 if (ao) {
1143 CHKERR AOApplicationToPetsc(ao, rind.size(), &*rind.begin());
1144 CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
1145 }
1146 auto it =
1147 set_block(add_bmc,
bit.rowField,
bit.colField,
bit.rowSide,
1148 bit.colSide,
bit.rowType,
bit.colType,
m, rind, cind);
1149 }
1150 }
1151
1152 for (auto &bit_col : bmc) {
1153 if (bit_col.setAtElement && bit_col.rowField == field &&
1154 bit_col.colField != field) {
1157 invMat.resize(inv.size1(), cm.size2(),
false);
1158 noalias(
invMat) = prod(inv, cm);
1159 if (ao)
1160 CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
1161 for (auto &bit_row : bmc) {
1162 if (bit_row.setAtElement && bit_row.rowField != field &&
1163 bit_row.colField == field) {
1166 K.resize(rind.size(), cind.size(),
false);
1167 noalias(
K) = prod(rm,
invMat);
1168 if (ao)
1169 CHKERR AOApplicationToPetsc(ao, rind.size(), &*rind.begin());
1170 auto it =
1171 add_block(add_bmc, bit_row.rowField, bit_col.colField,
1172 bit_row.rowSide, bit_col.colSide, bit_row.rowType,
1173 bit_col.colType,
K, rind, cind);
1174 }
1175 }
1176 }
1177 }
1178
1180 };
1181
1183 bool debug =
false) {
1186 for (
auto &
bit : add_bmc) {
1187 if (
bit.setAtElement) {
1188 std::cerr <<
"assemble vol: " <<
bit.rowField <<
" "
1189 <<
bit.colField << endl;
1190
1191 }
1192 }
1193 std::cerr << std::endl;
1194 }
1195 for (
auto &
bit : add_bmc) {
1196 if (
bit.setAtElement)
1198 }
1200 };
1201
1204 const std::string field,
1208
1209 for (
auto &
bit : add_bmc) {
1210 bit.unSetAtElement();
1212 }
1213
1214 for (
auto &
bit : bmc) {
1215 if (
bit.setAtElement) {
1216 if (
bit.rowField != field ||
bit.colField != field)
1217 auto it = set_block(add_bmc,
bit.rowField,
bit.colField,
1220 }
1221 }
1222
1224 bmc.get<1>().find(boost::make_tuple(field, field, row_type, row_type));
1225 if (
bit->setAtElement &&
bit != bmc.get<1>().end()) {
1226 auto it = set_block(add_bmc,
bit->rowField,
bit->colField,
1227 bit->rowSide,
bit->colSide,
bit->rowType,
1230 m +=
eps * diag_mat;
1231 } else {
1232 auto row_it = bmc.get<3>().lower_bound(field);
1233 for (; row_it != bmc.get<3>().upper_bound(field); ++row_it) {
1234 if (row_it->setAtElement) {
1235 auto it = set_block(add_bmc, field, field, 0, 0, row_type, row_type,
1236 diag_mat, row_it->rowInd, row_it->rowInd);
1239 break;
1240 }
1241 }
1242 if (row_it == bmc.get<3>().end())
1244 "row field not found %s", field.c_str());
1245 }
1246
1248 };
1249
1250 const bool debug =
false;
1252
1253 CHKERR invert_symm_schur(mat_cont,
"EP",
1255 CHKERR create_block_schur(mat_cont,
1258
1260 constexpr double eps = 1e-8;
1262 mass_mat_scalar,
eps);
1263
1264
1265
1266
1267
1268
1269
1271
1272 if(STauTau){
1278 }
1279
1281
1282 }
1283 }
1284
1286}
#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_OPERATION_UNSUCCESSFUL
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
FTensor::Index< 'm', SPACE_DIM > m
static __CLPK_integer lapack_dsysv(char uplo, __CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb, __CLPK_doublereal *work, __CLPK_integer lwork)
UBlasVector< int > VectorInt
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
multi_index_container< BlockMatData, indexed_by< ordered_unique< composite_key< BlockMatData, member< BlockMatData, std::string, &BlockMatData::rowField >, member< BlockMatData, std::string, &BlockMatData::colField >, member< BlockMatData, EntityType, &BlockMatData::rowType >, member< BlockMatData, EntityType, &BlockMatData::colType >, member< BlockMatData, int, &BlockMatData::rowSide >, member< BlockMatData, int, &BlockMatData::colSide > > >, ordered_non_unique< composite_key< BlockMatData, member< BlockMatData, std::string, &BlockMatData::rowField >, member< BlockMatData, std::string, &BlockMatData::colField >, member< BlockMatData, EntityType, &BlockMatData::rowType >, member< BlockMatData, EntityType, &BlockMatData::colType > > >, ordered_non_unique< composite_key< BlockMatData, member< BlockMatData, std::string, &BlockMatData::rowField >, member< BlockMatData, std::string, &BlockMatData::colField > > >, ordered_non_unique< member< BlockMatData, std::string, &BlockMatData::rowField > >, ordered_non_unique< member< BlockMatData, std::string, &BlockMatData::colField > > > > BlockMatContainer