v0.14.0
Loading...
Searching...
No Matches
Public Member Functions | Private Attributes | List of all members
OpLhsSkeleton Struct Reference

Operator the left hand side matrix. More...

Inheritance diagram for OpLhsSkeleton:
[legend]
Collaboration diagram for OpLhsSkeleton:
[legend]

Public Member Functions

 OpLhsSkeleton (boost::shared_ptr< FaceSideEle > side_fe_ptr, boost::shared_ptr< MatrixDouble > d_mat_fluid_ptr, boost::shared_ptr< MatrixDouble > d_mat_solid_ptr, boost::shared_ptr< MatrixDouble > d_mat_visc_fluid_ptr, boost::shared_ptr< MatrixDouble > d_mat_visc_solid_ptr, boost::shared_ptr< Range > r_visc_fluid_ptr, boost::shared_ptr< Range > r_visc_solid_ptr, SmartPetscObj< Mat > K, SmartPetscObj< Mat > C, SmartPetscObj< Mat > M)
 Construct a new OpH1LhsSkeleton. More...
 
MoFEMErrorCode doWork (int side, EntityType type, EntitiesFieldData::EntData &data)
 

Private Attributes

boost::shared_ptr< FaceSideElesideFEPtr
 pointer to element to get data on edge/face sides More...
 
MatrixDouble locMatK
 
MatrixDouble locMatC
 
MatrixDouble locMatM
 
boost::shared_ptr< MatrixDouble > dMatFluidPtr
 
boost::shared_ptr< MatrixDouble > dMatSolidPtr
 
boost::shared_ptr< MatrixDouble > dViscFluidPtr
 
boost::shared_ptr< MatrixDouble > dViscSolidPtr
 
boost::shared_ptr< RangerFluidPtr
 
boost::shared_ptr< RangerSolidPtr
 
SmartPetscObj< Mat > K
 
SmartPetscObj< Mat > C
 
SmartPetscObj< Mat > M
 

Detailed Description

Operator the left hand side matrix.

Examples
level_set.cpp.

Definition at line 110 of file fluid_structure_eigenproblem.cpp.

Constructor & Destructor Documentation

◆ OpLhsSkeleton()

OpLhsSkeleton::OpLhsSkeleton ( boost::shared_ptr< FaceSideEle side_fe_ptr,
boost::shared_ptr< MatrixDouble >  mat_d_fluid_ptr,
boost::shared_ptr< MatrixDouble >  mat_d_solid_ptr,
boost::shared_ptr< MatrixDouble >  mat_d_visc_fluid_ptr,
boost::shared_ptr< MatrixDouble >  mat_d_visc_solid_ptr,
boost::shared_ptr< Range r_fluid_ptr,
boost::shared_ptr< Range r_solid_ptr,
SmartPetscObj< Mat >  K,
SmartPetscObj< Mat >  C,
SmartPetscObj< Mat >  M 
)

Construct a new OpH1LhsSkeleton.

Parameters
side_fe_ptrpointer to FE to evaluate side elements
Parameters
side_fe_ptrpointer to FE to evaluate side elements

Definition at line 1042 of file fluid_structure_eigenproblem.cpp.

1050 : SkeletonEleOp(NOSPACE, SkeletonEleOp::OPSPACE), sideFEPtr(side_fe_ptr),
1051 dMatFluidPtr(mat_d_fluid_ptr), dMatSolidPtr(mat_d_solid_ptr),
1052 dViscFluidPtr(mat_d_visc_fluid_ptr), dViscSolidPtr(mat_d_visc_solid_ptr),
1053 rFluidPtr(r_fluid_ptr), rSolidPtr(r_solid_ptr), K(K), C(C), M(M) {}
@ NOSPACE
Definition: definitions.h:83
ElementsAndOps< SPACE_DIM >::SkeletonEleOp SkeletonEleOp
boost::shared_ptr< MatrixDouble > dMatSolidPtr
boost::shared_ptr< Range > rSolidPtr
boost::shared_ptr< MatrixDouble > dMatFluidPtr
boost::shared_ptr< Range > rFluidPtr
boost::shared_ptr< MatrixDouble > dViscSolidPtr
boost::shared_ptr< FaceSideEle > sideFEPtr
pointer to element to get data on edge/face sides
boost::shared_ptr< MatrixDouble > dViscFluidPtr

Member Function Documentation

◆ doWork()

MoFEMErrorCode OpLhsSkeleton::doWork ( int  side,
EntityType  type,
EntitiesFieldData::EntData data 
)

Definition at line 1055 of file fluid_structure_eigenproblem.cpp.

1056 {
1058
1059 // get normal of the face or edge
1060 auto t_normal = getFTensor1Normal();
1061 t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
1062
1063 // Collect data from side domian elements
1064 CHKERR loopSide("dFE", sideFEPtr.get(), SPACE_DIM);
1065 const auto in_the_loop =
1066 sideFEPtr->nInTheLoop; // return number of elements on the side
1067
1068 std::array<boost::shared_ptr<MatrixDouble>, 2> d_ptr;
1069 std::array<boost::shared_ptr<MatrixDouble>, 2> d_vis_ptr;
1070
1071 std::array<std::vector<VectorInt>, 2> indices_side_map;
1072 std::array<std::vector<MatrixDouble>, 2> base_side_map;
1073 std::array<std::vector<MatrixDouble>, 2> diff_base_side_map;
1074
1075 // MOFEM_LOG("SELF", Sev::inform)
1076 // << "AAA " << sideHandle[0] << " " << sideHandle[1];
1077
1078 // auto save_range = [&](const std::string name, const Range r) {
1079 // MoFEMFunctionBegin;
1080 // EntityHandle out_meshset;
1081 // CHKERR getPtrFE()->mField.get_moab().create_meshset(MESHSET_SET,
1082 // out_meshset);
1083 // CHKERR getPtrFE()->mField.get_moab().add_entities(out_meshset, r);
1084 // CHKERR getPtrFE()->mField.get_moab().write_file(name.c_str(), "VTK", "",
1085 // &out_meshset, 1);
1086 // CHKERR getPtrFE()->mField.get_moab().delete_entities(&out_meshset, 1);
1087 // MoFEMFunctionReturn(0);
1088 // };
1089
1090 // Tag th_n;
1091 // double def[] = {0, 0, 0};
1092 // CHKERR getPtrFE()->mField.get_moab().tag_get_handle(
1093 // "NORMAL", 3, MB_TYPE_DOUBLE, th_n, MB_TAG_CREAT | MB_TAG_SPARSE, def);
1094
1095 auto set_row_col = [&](auto a, auto b) {
1096 d_ptr[a] = dMatSolidPtr;
1097 d_ptr[b] = dMatFluidPtr;
1098 d_vis_ptr[a] = dViscSolidPtr;
1099 d_vis_ptr[b] = dViscFluidPtr;
1100 indices_side_map[a].swap(indicesUsSideMap[a]);
1101 indices_side_map[b].swap(indicesUfSideMap[b]);
1102 base_side_map[a].swap(baseUsSideMap[a]);
1103 base_side_map[b].swap(baseUfSideMap[b]);
1104 diff_base_side_map[a].swap(diffUsBaseSideMap[a]);
1105 diff_base_side_map[b].swap(diffUfBaseSideMap[b]);
1106 };
1107
1108 if (rSolidPtr->find(sideHandle[0]) != rSolidPtr->end()) {
1109
1110 if (rFluidPtr->find(sideHandle[1]) == rFluidPtr->end()) {
1111 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1112 "Other side should be fluid");
1113 }
1114
1115 set_row_col(0, 1);
1116
1117 if (senseMap[0] < 0)
1118 t_normal(i) *= -1;
1119
1120 } else if (rFluidPtr->find(sideHandle[0]) != rFluidPtr->end()) {
1121
1122 if (rSolidPtr->find(sideHandle[1]) == rSolidPtr->end()) {
1123 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1124 "Other side should be solid");
1125 }
1126
1127 set_row_col(1, 0);
1128
1129 if (senseMap[1] < 0)
1130 t_normal(i) *= -1;
1131
1132 } else {
1133 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1134 "Side is not fluid or solid");
1135 }
1136
1137 // auto fe_h = getFEEntityHandle();
1138 // CHKERR getPtrFE()->mField.get_moab().tag_set_data(th_n, &fe_h, 1,
1139 // &t_normal(0));
1140
1141 // calculate penalty
1142 const double s = getMeasure() / (areaMap[0] + areaMap[1]);
1143 const double p = pow(penalty, order) * s;
1144
1145 // get number of integration points on face
1146 const size_t nb_integration_pts = getGaussPts().size2();
1147
1148 // when boundary edge/face is evaluated, and 2 when is skeleton edge/face.
1149 const double beta = 1. / (in_the_loop + 1);
1150
1151 auto integrate = [&](auto sense_row, auto &row_ind, auto &row, auto &row_diff,
1152 auto sense_col, auto &col_ind, auto &col, auto &col_diff,
1153 auto &d_row, auto &d_col, auto &d_vis_row,
1154 auto &d_vis_col) {
1156
1157 // number of dofs, for homogenous approximation this value is
1158 // constant.
1159 const auto nb_rows = row_ind.size();
1160 const auto nb_cols = col_ind.size();
1161
1162 const auto nb_row_base_functions = row.size2();
1163
1164 auto t_D_row = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(d_row);
1165 auto t_D_col = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(d_col);
1166
1167 auto t_vis_D_row =
1168 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(d_vis_row);
1169 auto t_vis_D_col =
1170 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(d_vis_col);
1171
1172 if (nb_cols) {
1173
1174 // resize local element matrix
1175 locMatK.resize(nb_rows, nb_cols, false);
1176 locMatK.clear();
1177 locMatC.resize(nb_rows, nb_cols, false);
1178 locMatC.clear();
1179 locMatM.resize(nb_rows, nb_cols, false);
1180 locMatM.clear();
1181
1182 // get base functions, and integration weights
1183 auto t_row_base = get_ntensor(row);
1184 auto t_diff_row_base = get_diff_ntensor(row_diff);
1185
1186 constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
1188 t_P(i, j) = t_normal(i) * t_normal(j);
1189
1190 auto t_w = getFTensor0IntegrationWeight();
1191
1192 // iterate integration points on face/edge
1193 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1194
1195 // t_w is integration weight, and measure is element area, or
1196 // volume, depending if problem is in 2d/3d.
1197 const double alpha = getMeasure() * t_w;
1198
1199 // iterate rows
1200 size_t rr = 0;
1201 for (; rr != nb_rows / SPACE_DIM; ++rr) {
1202
1203 auto get_vn_plus = [&](auto &t_D) {
1205 t_vn_plus(i, j, k) =
1206 (beta * phi / p) *
1207 (t_kd(i, m) * t_D(m, j, k, l) * t_diff_row_base(l));
1208 return t_vn_plus;
1209 };
1210
1211 auto get_vn = [&]() {
1213 t_vn(i, j, k) = t_normal(j) * (t_kd(i, k) * t_row_base * sense_row);
1214 return t_vn;
1215 };
1216
1217 auto t_vn_plus = get_vn_plus(t_D_row);
1218 auto t_vn = get_vn();
1219 auto t_vn_lambda_plus = get_vn_plus(t_vis_D_row);
1220
1221 // get base functions on columns
1222 auto t_col_base = get_ntensor(col, gg, 0);
1223 auto t_diff_col_base = get_diff_ntensor(col_diff, gg, 0);
1224
1225 auto t_mat_K = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
1226 locMatK, SPACE_DIM * rr);
1227 auto t_mat_C = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
1228 locMatC, SPACE_DIM * rr);
1229 auto t_mat_M = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
1230 locMatM, SPACE_DIM * rr);
1231
1232 // iterate columns
1233 for (size_t cc = 0; cc != nb_cols / SPACE_DIM; ++cc) {
1234
1235 auto get_un_plus = [&](auto &t_D) {
1237 t_un_plus;
1238 t_un_plus(i, j, k) =
1239 beta * ((t_P(i, m) * t_D(m, j, k, l)) * t_diff_col_base(l));
1240 return t_un_plus;
1241 };
1242
1243 auto get_un = [&]() {
1245 t_un(i, j, k) =
1246 (t_normal(j) * (t_P(i, k) * t_col_base * sense_col));
1247 return t_un;
1248 };
1249
1250 auto t_un_plus = get_un_plus(t_D_col);
1251 auto t_un = get_un();
1252 auto t_un_lambda_plus = get_un_plus(t_vis_D_col);
1253
1254 // K
1255 t_mat_K(i, j) -=
1256 alpha * (t_vn(m, n, i) - t_vn_plus(m, n, i)) *
1257 ((-p) * (t_un(m, n, j) - (t_un_plus(m, n, j) / p)));
1258 t_mat_K(i, j) -= alpha * (t_vn_plus(m, n, i) * t_un_plus(m, n, j));
1259
1260 // C1
1261 t_mat_C(i, j) -=
1262 alpha * (-t_vn_lambda_plus(m, n, i)) *
1263 ((-p) * (t_un(m, n, j) - (t_un_plus(m, n, j) / p)));
1264 t_mat_C(i, j) -=
1265 alpha * (t_vn_lambda_plus(m, n, i) * t_un_plus(m, n, j));
1266
1267 // C2
1268 t_mat_C(i, j) -= alpha * (t_vn(m, n, i) - t_vn_plus(m, n, i)) *
1269 ((-p) * (t_un_lambda_plus(m, n, j) / (-p)));
1270 t_mat_C(i, j) -=
1271 alpha * (t_vn_plus(m, n, i) * t_un_lambda_plus(m, n, j));
1272
1273 // M
1274 t_mat_M(i, j) -= alpha * (-t_vn_lambda_plus(m, n, i)) *
1275 ((-p) * (t_un_lambda_plus(m, n, j) / (-p)));
1276 t_mat_M(i, j) -=
1277 alpha * (t_vn_lambda_plus(m, n, i) * t_un_lambda_plus(m, n, j));
1278
1279 // move to next column base and element of matrix
1280 ++t_col_base;
1281 ++t_diff_col_base;
1282 ++t_mat_K;
1283 ++t_mat_C;
1284 ++t_mat_M;
1285 }
1286
1287 // move to next row base
1288 ++t_row_base;
1289 ++t_diff_row_base;
1290 }
1291
1292 for (; rr < nb_row_base_functions; ++rr) {
1293 ++t_row_base;
1294 ++t_diff_row_base;
1295 }
1296
1297 ++t_w;
1298 }
1299
1300 // assemble system
1301
1302 CHKERR ::MatSetValues(K, nb_rows, &*row_ind.begin(), col_ind.size(),
1303 &*col_ind.begin(), &*locMatK.data().begin(),
1304 ADD_VALUES);
1305 CHKERR ::MatSetValues(C, nb_rows, &*row_ind.begin(), col_ind.size(),
1306 &*col_ind.begin(), &*locMatC.data().begin(),
1307 ADD_VALUES);
1308 CHKERR ::MatSetValues(M, nb_rows, &*row_ind.begin(), col_ind.size(),
1309 &*col_ind.begin(), &*locMatM.data().begin(),
1310 ADD_VALUES);
1311 }
1312
1314 };
1315
1316 // iterate the sides rows
1317 for (auto s0 : {LEFT_SIDE, RIGHT_SIDE}) {
1318 const auto sense_row = senseMap[s0];
1319
1320 for (auto x0 = 0; x0 != indices_side_map[s0].size(); ++x0) {
1321
1322 for (auto s1 : {LEFT_SIDE, RIGHT_SIDE}) {
1323 const auto sense_col = senseMap[s1];
1324
1325 for (auto x1 = 0; x1 != indices_side_map[s1].size(); ++x1) {
1326
1327 CHKERR integrate(sense_row, indices_side_map[s0][x0],
1328 base_side_map[s0][x0], diff_base_side_map[s0][x0],
1329
1330 sense_col, indices_side_map[s1][x1],
1331 base_side_map[s1][x1], diff_base_side_map[s1][x1],
1332
1333 *d_ptr[s0], *d_ptr[s1], *d_vis_ptr[s0],
1334 *d_vis_ptr[s1]
1335
1336 );
1337 }
1338 }
1339 }
1340 }
1341
1343}
static Index< 'p', 3 > p
constexpr double a
Kronecker Delta class.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
auto get_diff_ntensor(T &base_mat)
std::array< std::vector< MatrixDouble >, 2 > baseUfSideMap
FTensor::Index< 'j', SPACE_DIM > j
FTensor::Index< 'k', SPACE_DIM > k
std::array< std::vector< MatrixDouble >, 2 > baseUsSideMap
FTensor::Index< 'i', SPACE_DIM > i
auto get_ntensor(T &base_mat)
constexpr int SPACE_DIM
FTensor::Index< 'l', SPACE_DIM > l
std::array< double, 2 > areaMap
std::array< int, 2 > senseMap
std::array< std::vector< VectorInt >, 2 > indicesUfSideMap
FTensor::Index< 'n', SPACE_DIM > n
std::array< EntityHandle, 2 > sideHandle
std::array< std::vector< MatrixDouble >, 2 > diffUfBaseSideMap
static double penalty
FTensor::Index< 'm', SPACE_DIM > m
std::array< std::vector< MatrixDouble >, 2 > diffUsBaseSideMap
static double phi
std::array< std::vector< VectorInt >, 2 > indicesUsSideMap
constexpr auto t_kd

Member Data Documentation

◆ C

SmartPetscObj<Mat> OpLhsSkeleton::C
private

Definition at line 142 of file fluid_structure_eigenproblem.cpp.

◆ dMatFluidPtr

boost::shared_ptr<MatrixDouble> OpLhsSkeleton::dMatFluidPtr
private

Definition at line 135 of file fluid_structure_eigenproblem.cpp.

◆ dMatSolidPtr

boost::shared_ptr<MatrixDouble> OpLhsSkeleton::dMatSolidPtr
private

Definition at line 136 of file fluid_structure_eigenproblem.cpp.

◆ dViscFluidPtr

boost::shared_ptr<MatrixDouble> OpLhsSkeleton::dViscFluidPtr
private

Definition at line 137 of file fluid_structure_eigenproblem.cpp.

◆ dViscSolidPtr

boost::shared_ptr<MatrixDouble> OpLhsSkeleton::dViscSolidPtr
private

Definition at line 138 of file fluid_structure_eigenproblem.cpp.

◆ K

SmartPetscObj<Mat> OpLhsSkeleton::K
private

Definition at line 141 of file fluid_structure_eigenproblem.cpp.

◆ locMatC

MatrixDouble OpLhsSkeleton::locMatC
private

Definition at line 133 of file fluid_structure_eigenproblem.cpp.

◆ locMatK

MatrixDouble OpLhsSkeleton::locMatK
private

Definition at line 132 of file fluid_structure_eigenproblem.cpp.

◆ locMatM

MatrixDouble OpLhsSkeleton::locMatM
private

Definition at line 134 of file fluid_structure_eigenproblem.cpp.

◆ M

SmartPetscObj<Mat> OpLhsSkeleton::M
private

Definition at line 143 of file fluid_structure_eigenproblem.cpp.

◆ rFluidPtr

boost::shared_ptr<Range> OpLhsSkeleton::rFluidPtr
private

Definition at line 139 of file fluid_structure_eigenproblem.cpp.

◆ rSolidPtr

boost::shared_ptr<Range> OpLhsSkeleton::rSolidPtr
private

Definition at line 140 of file fluid_structure_eigenproblem.cpp.

◆ sideFEPtr

boost::shared_ptr<FaceSideEle> OpLhsSkeleton::sideFEPtr
private

pointer to element to get data on edge/face sides

Definition at line 131 of file fluid_structure_eigenproblem.cpp.


The documentation for this struct was generated from the following file: