v0.15.0
Loading...
Searching...
No Matches
Public Member Functions | Public Attributes | List of all members
ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumJacobian Struct Reference
User modules » Mass Element

#include "users_modules/basic_finite_elements/src/ConvectiveMassElement.hpp"

Inheritance diagram for ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumJacobian:
[legend]
Collaboration diagram for ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumJacobian:
[legend]

Public Member Functions

 OpEshelbyDynamicMaterialMomentumJacobian (const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian=true)
 
MoFEMErrorCode doWork (int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
 

Public Attributes

BlockData & dAta
 
CommonData & commonData
 
int tAg
 
bool jAcobian
 
bool fieldDisp
 
VectorBoundedArray< adouble, 3 > a
 
VectorBoundedArray< adouble, 3 > v
 
VectorBoundedArray< adouble, 3 > a_T
 
MatrixBoundedArray< adouble, 9 > g
 
MatrixBoundedArray< adouble, 9 > H
 
MatrixBoundedArray< adouble, 9 > invH
 
MatrixBoundedArray< adouble, 9 > h
 
MatrixBoundedArray< adouble, 9 > F
 
MatrixBoundedArray< adouble, 9 > G
 
VectorDouble active
 

Detailed Description

Definition at line 347 of file ConvectiveMassElement.hpp.

Constructor & Destructor Documentation

◆ OpEshelbyDynamicMaterialMomentumJacobian()

ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumJacobian::OpEshelbyDynamicMaterialMomentumJacobian ( const std::string field_name,
BlockData & data,
CommonData & common_data,
int tag,
bool jacobian = true )

Definition at line 1234 of file ConvectiveMassElement.cpp.

1239 : VolumeElementForcesAndSourcesCore::UserDataOperator(
1240 field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
1241 dAta(data), commonData(common_data), tAg(tag), jAcobian(jacobian),
1242 fieldDisp(false) {}
field_name
constexpr auto field_name
Definition poisson_2d_homogeneous.cpp:13

Member Function Documentation

◆ doWork()

MoFEMErrorCode ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumJacobian::doWork ( int row_side,
EntityType row_type,
EntitiesFieldData::EntData & row_data )

Definition at line 1245 of file ConvectiveMassElement.cpp.

1246 {
1247 MoFEMFunctionBeginHot;
1248
1249 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1250 dAta.tEts.end()) {
1251 MoFEMFunctionReturnHot(0);
1252 }
1253
1254 // do it only once, no need to repeat this for edges,faces or tets
1255 if (row_type != MBVERTEX)
1256 MoFEMFunctionReturnHot(0);
1257
1258 int nb_dofs = row_data.getIndices().size();
1259 if (nb_dofs == 0)
1260 MoFEMFunctionReturnHot(0);
1261
1262 try {
1263
1264 a.resize(3);
1265 v.resize(3);
1266 g.resize(3, 3);
1267 G.resize(3, 3);
1268 h.resize(3, 3);
1269 F.resize(3, 3);
1270 H.resize(3, 3);
1271 H.clear();
1272 invH.resize(3, 3);
1273 invH.clear();
1274 for (int dd = 0; dd < 3; dd++) {
1275 H(dd, dd) = 1;
1276 invH(dd, dd) = 1;
1277 }
1278
1279 int nb_gauss_pts = row_data.getN().size1();
1280 commonData.valT.resize(nb_gauss_pts);
1281 commonData.jacTRowPtr.resize(nb_gauss_pts);
1282 commonData.jacT.resize(nb_gauss_pts);
1283
1284 int nb_active_vars = 0;
1285 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1286
1287 if (gg == 0) {
1288
1289 trace_on(tAg);
1290
1291 for (int nn1 = 0; nn1 < 3; nn1++) { // 0
1292 a[nn1] <<=
1293 commonData.dataAtGaussPts["DOT_" + commonData.spatialVelocities]
1294 [gg][nn1];
1295 nb_active_vars++;
1296 }
1297
1298 for (int nn1 = 0; nn1 < 3; nn1++) { // 3
1299 v[nn1] <<=
1300 commonData.dataAtGaussPts[commonData.spatialVelocities][gg][nn1];
1301 nb_active_vars++;
1302 }
1303 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3
1304 for (int nn2 = 0; nn2 < 3; nn2++) {
1305 g(nn1, nn2) <<=
1306 commonData.gradAtGaussPts[commonData.spatialVelocities][gg](
1307 nn1, nn2);
1308 nb_active_vars++;
1309 }
1310 }
1311 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9
1312 for (int nn2 = 0; nn2 < 3; nn2++) {
1313 h(nn1, nn2) <<=
1314 commonData.gradAtGaussPts[commonData.spatialPositions][gg](nn1,
1315 nn2);
1316 nb_active_vars++;
1317 if (fieldDisp) {
1318 if (nn1 == nn2) {
1319 h(nn1, nn2) += 1;
1320 }
1321 }
1322 }
1323 }
1324 if (commonData.gradAtGaussPts[commonData.meshPositions].size() > 0) {
1325 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9+9
1326 for (int nn2 = 0; nn2 < 3; nn2++) {
1327 H(nn1, nn2) <<=
1328 commonData.gradAtGaussPts[commonData.meshPositions][gg](nn1,
1329 nn2);
1330 nb_active_vars++;
1331 }
1332 }
1333 }
1334 adouble detH;
1335 detH = 1;
1336 if (commonData.gradAtGaussPts[commonData.meshPositions].size() > 0) {
1337 detH = determinantTensor3by3(H);
1338 }
1339 CHKERR invertTensor3by3(H, detH, invH);
1340
1341 FTensor::Index<'i', 3> i;
1342 FTensor::Index<'j', 3> j;
1343 FTensor::Index<'k', 3> k;
1344
1345 a_T.resize(3);
1346
1347 auto t_h = getFTensor2FromArray3by3(h, FTensor::Number<0>(), 0);
1348 auto t_invH = getFTensor2FromArray3by3(invH, FTensor::Number<0>(), 0);
1349 auto t_F = getFTensor2FromArray3by3(F, FTensor::Number<0>(), 0);
1350 auto t_g = getFTensor2FromArray3by3(g, FTensor::Number<0>(), 0);
1351 auto t_G = getFTensor2FromArray3by3(G, FTensor::Number<0>(), 0);
1352
1353 auto t_a = FTensor::Tensor1<adouble *, 3>{&a[0], &a[1], &a[2]};
1354 auto t_v = FTensor::Tensor1<adouble *, 3>{&v[0], &v[1], &v[2]};
1355 auto t_a_T = FTensor::Tensor1<adouble *, 3>{&a_T[0], &a_T[1], &a_T[2]};
1356
1357 t_F(i, j) = t_h(i, k) * t_invH(k, j);
1358 t_G(i, j) = t_g(i, k) * t_invH(k, j);
1359 t_a_T(i) = t_F(k, i) * t_a(k) + t_G(k, i) * t_v(k);
1360 const auto rho0 = dAta.rho0;
1361 t_a_T(i) = -rho0 * detH;
1362
1363 commonData.valT[gg].resize(3);
1364 for (int nn = 0; nn < 3; nn++) {
1365 a_T[nn] >>= (commonData.valT[gg])[nn];
1366 }
1367 trace_off();
1368 }
1369
1370 active.resize(nb_active_vars);
1371 int aa = 0;
1372 for (int nn1 = 0; nn1 < 3; nn1++) { // 0
1373 active[aa++] =
1374 commonData
1375 .dataAtGaussPts["DOT_" + commonData.spatialVelocities][gg][nn1];
1376 }
1377
1378 for (int nn1 = 0; nn1 < 3; nn1++) { // 3
1379 active[aa++] =
1380 commonData.dataAtGaussPts[commonData.spatialVelocities][gg][nn1];
1381 }
1382 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3
1383 for (int nn2 = 0; nn2 < 3; nn2++) {
1384 active[aa++] =
1385 commonData.gradAtGaussPts[commonData.spatialVelocities][gg](nn1,
1386 nn2);
1387 }
1388 }
1389 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9
1390 for (int nn2 = 0; nn2 < 3; nn2++) {
1391 if (fieldDisp && nn1 == nn2) {
1392 active[aa++] =
1393 commonData.gradAtGaussPts[commonData.spatialPositions][gg](
1394 nn1, nn2) +
1395 1;
1396 } else {
1397 active[aa++] =
1398 commonData.gradAtGaussPts[commonData.spatialPositions][gg](nn1,
1399 nn2);
1400 }
1401 }
1402 }
1403 if (commonData.gradAtGaussPts[commonData.meshPositions].size() > 0) {
1404 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9+9
1405 for (int nn2 = 0; nn2 < 3; nn2++) {
1406 active[aa++] =
1407 commonData.gradAtGaussPts[commonData.meshPositions][gg](nn1,
1408 nn2);
1409 }
1410 }
1411 }
1412
1413 if (!jAcobian) {
1414 VectorDouble &res = commonData.valT[gg];
1415 if (gg > 0) {
1416 res.resize(3);
1417 int r;
1418 r = ::function(tAg, 3, nb_active_vars, &active[0], &res[0]);
1419 if (r != 3) { // function is locally analytic
1420 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1421 "ADOL-C function evaluation with error r = %d", r);
1422 }
1423 }
1424 double val = getVolume() * getGaussPts()(3, gg);
1425 res *= val;
1426 } else {
1427 commonData.jacTRowPtr[gg].resize(3);
1428 commonData.jacT[gg].resize(3, nb_active_vars);
1429 for (int nn1 = 0; nn1 < 3; nn1++) {
1430 (commonData.jacTRowPtr[gg])[nn1] = &(commonData.jacT[gg](nn1, 0));
1431 }
1432 int r;
1433 r = jacobian(tAg, 3, nb_active_vars, &active[0],
1434 &(commonData.jacTRowPtr[gg])[0]);
1435 if (r != 3) {
1436 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1437 "ADOL-C function evaluation with error");
1438 }
1439 double val = getVolume() * getGaussPts()(3, gg);
1440 commonData.jacT[gg] *= val;
1441 }
1442 }
1443
1444 } catch (const std::exception &ex) {
1445 std::ostringstream ss;
1446 ss << "throw in method: " << ex.what() << std::endl;
1447 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "%s", ss.str().c_str());
1448 }
1449
1450 MoFEMFunctionReturnHot(0);
1451}
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:456
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:19
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:20
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition ddTensor0.hpp:33
MoFEM::invertTensor3by3
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition Templates.hpp:1622
MoFEM::getFTensor2FromArray3by3
auto getFTensor2FromArray3by3(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr, const size_t cc=0)
Definition Templates.hpp:1376
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition Templates.hpp:1574
sdf.r
int r
Definition sdf.py:8
ConvectiveMassElement::CommonData::jacTRowPtr
std::vector< std::vector< double * > > jacTRowPtr
Definition ConvectiveMassElement.hpp:145
ConvectiveMassElement::CommonData::dataAtGaussPts
std::map< std::string, std::vector< VectorDouble > > dataAtGaussPts
Definition ConvectiveMassElement.hpp:133
ConvectiveMassElement::CommonData::gradAtGaussPts
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
Definition ConvectiveMassElement.hpp:134
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition EntitiesFieldData.hpp:1318
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition EntitiesFieldData.hpp:1214

Member Data Documentation

◆ a

VectorBoundedArray<adouble, 3> ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumJacobian::a

Definition at line 362 of file ConvectiveMassElement.hpp.

◆ a_T

VectorBoundedArray<adouble, 3> ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumJacobian::a_T

Definition at line 362 of file ConvectiveMassElement.hpp.

◆ active

VectorDouble ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumJacobian::active

Definition at line 364 of file ConvectiveMassElement.hpp.

◆ commonData

CommonData& ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumJacobian::commonData

Definition at line 352 of file ConvectiveMassElement.hpp.

◆ dAta

BlockData& ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumJacobian::dAta

Definition at line 351 of file ConvectiveMassElement.hpp.

◆ F

MatrixBoundedArray<adouble, 9> ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumJacobian::F

Definition at line 363 of file ConvectiveMassElement.hpp.

◆ fieldDisp

bool ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumJacobian::fieldDisp

Definition at line 355 of file ConvectiveMassElement.hpp.

◆ G

MatrixBoundedArray<adouble, 9> ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumJacobian::G

Definition at line 363 of file ConvectiveMassElement.hpp.

◆ g

MatrixBoundedArray<adouble, 9> ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumJacobian::g

Definition at line 363 of file ConvectiveMassElement.hpp.

◆ H

MatrixBoundedArray<adouble, 9> ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumJacobian::H

Definition at line 363 of file ConvectiveMassElement.hpp.

◆ h

MatrixBoundedArray<adouble, 9> ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumJacobian::h

Definition at line 363 of file ConvectiveMassElement.hpp.

◆ invH

MatrixBoundedArray<adouble, 9> ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumJacobian::invH

Definition at line 363 of file ConvectiveMassElement.hpp.

◆ jAcobian

bool ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumJacobian::jAcobian

Definition at line 354 of file ConvectiveMassElement.hpp.

◆ tAg

int ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumJacobian::tAg

Definition at line 353 of file ConvectiveMassElement.hpp.

◆ v

VectorBoundedArray<adouble, 3> ConvectiveMassElement::OpEshelbyDynamicMaterialMomentumJacobian::v

Definition at line 362 of file ConvectiveMassElement.hpp.

The documentation for this struct was generated from the following files:
Generated by Doxygen 1.12.0 and hosted at University of Glasgow