Operator for linear form, usually to calculate values on right hand side.
292 {
294
295 if (row_type != MBVERTEX) {
297 }
298
300
301 auto tensor_to_tensor = [](const auto &t1, auto &t2) {
302 t2(0, 0) = t1(0, 0);
303 t2(1, 1) = t1(1, 1);
304 t2(2, 2) = t1(2, 2);
305 t2(0, 1) = t2(1, 0) = t1(1, 0);
306 t2(0, 2) = t2(2, 0) = t1(2, 0);
307 t2(1, 2) = t2(2, 1) = t1(2, 1);
308 };
309
310 auto tensor_to_vector = [](
const auto &
t,
auto &
v) {
317 };
318
319 auto get_tag_handle = [&](auto name, auto size) {
321 std::vector<double> def_vals(size, 0.0);
323 MB_TAG_CREAT | MB_TAG_SPARSE,
324 def_vals.data());
326 };
327
328 const int nb_integration_pts = row_data.getN().size1();
329
330 auto th_internal_stress =
331 get_tag_handle("INTERNAL_STRESS", 9 * nb_integration_pts);
332 auto th_actual_stress =
333 get_tag_handle("ACTUAL_STRESS", 9 * nb_integration_pts);
334 auto th_deviatoric_stress =
335 get_tag_handle("DEVIATORIC_STRESS", 9 * nb_integration_pts);
336 auto th_hydrostatic_stress =
337 get_tag_handle("HYDROSTATIC_STRESS", 9 * nb_integration_pts);
338
339 auto th_internal_stress_mean = get_tag_handle("MED_INTERNAL_STRESS", 9);
340 auto th_actual_stress_mean = get_tag_handle("MED_ACTUAL_STRESS", 9);
341
346
347 auto t_h = getFTensor2FromMat<3, 3>(*
dataAtPts->hMat);
348 auto t_H = getFTensor2FromMat<3, 3>(*
dataAtPts->HMat);
349
350 dataAtPts->stiffnessMat->resize(36, 1,
false);
354 const double young =
m.second.E;
355 const double poisson =
m.second.PoissonRatio;
356
357 const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
358
359 t_D(
i,
j,
k,
l) = 0.;
360 t_D(0, 0, 0, 0) = t_D(1, 1, 1, 1) = t_D(2, 2, 2, 2) = 1 - poisson;
361 t_D(0, 1, 0, 1) = t_D(0, 2, 0, 2) = t_D(1, 2, 1, 2) =
362 0.5 * (1 - 2 * poisson);
363 t_D(0, 0, 1, 1) = t_D(1, 1, 0, 0) = t_D(0, 0, 2, 2) = t_D(2, 2, 0, 0) =
364 t_D(1, 1, 2, 2) = t_D(2, 2, 1, 1) = poisson;
365 t_D(
i,
j,
k,
l) *= coefficient;
366
367 break;
368 }
369
370 double detH = 0.;
373
376
378
379 auto t_internal_stress =
380 getFTensor2FromMat<3, 3>(*
dataAtPts->internalStressMat);
381
382 dataAtPts->actualStressMat->resize(9, nb_integration_pts,
false);
383 auto t_actual_stress = getFTensor2FromMat<3, 3>(*
dataAtPts->actualStressMat);
384
385 dataAtPts->deviatoricStressMat->resize(9, nb_integration_pts,
false);
386 auto t_deviatoric_stress =
387 getFTensor2FromMat<3, 3>(*
dataAtPts->deviatoricStressMat);
388 dataAtPts->hydrostaticStressMat->resize(9, nb_integration_pts,
false);
389 auto t_hydrostatic_stress =
390 getFTensor2FromMat<3, 3>(*
dataAtPts->hydrostaticStressMat);
391
394
395 t_internal_stress_mean(
i,
j) = 0.;
396 t_actual_stress_mean(
i,
j) = 0.;
397
399
400 for (int gg = 0; gg != nb_integration_pts; ++gg) {
401
403 t_small_strain_symm(
i,
j) = (t_h(
i,
j) || t_h(
j,
i)) / 2.;
404 } else {
407 t_F(
i,
j) = t_h(
i,
k) * t_invH(
k,
j);
408 t_small_strain_symm(
i,
j) = (t_F(
i,
j) || t_F(
j,
i)) / 2.;
409 ++t_H;
410 }
411
413 for (auto ii : {0, 1, 2}) {
414 t_small_strain_symm(ii, ii) -= 1.;
415 }
416 }
417
418
419 t_stress_symm(
i,
j) = t_D(
i,
j,
k,
l) * t_small_strain_symm(
k,
l);
420 tensor_to_tensor(t_stress_symm, t_stress);
421 t_actual_stress(
i,
j) = t_stress(
i,
j) + t_internal_stress(
i,
j);
422
425
426 double hydrostatic_pressure =
427 (t_actual_stress(0, 0) + t_actual_stress(1, 1) +
428 t_actual_stress(2, 2)) /
429 3.;
430
431 t_hydrostatic_stress(
i,
j) = 0.;
432 for (auto ii : {0, 1, 2}) {
433 t_hydrostatic_stress(ii, ii) += hydrostatic_pressure;
434 }
435
436 t_deviatoric_stress(
i,
j) = t_actual_stress(
i,
j);
437 for (auto ii : {0, 1, 2}) {
438 t_deviatoric_stress(ii, ii) -= hydrostatic_pressure;
439 }
440
441 t_actual_stress_mean(
i,
j) += t_w * t_actual_stress(
i,
j);
442 t_internal_stress_mean(
i,
j) += t_w * t_internal_stress(
i,
j);
443
444 ++t_w;
445 ++t_h;
446 ++t_actual_stress;
447 ++t_internal_stress;
448 }
449
452 vec_actual_stress_mean.resize(9, false);
453 vec_actual_stress_mean.clear();
454
456 vec_internal_stress_mean.resize(9, false);
457 vec_internal_stress_mean.clear();
458
459 tensor_to_vector(t_actual_stress_mean, vec_actual_stress_mean);
460 tensor_to_vector(t_internal_stress_mean, vec_internal_stress_mean);
461
463 &*(vec_actual_stress_mean.data().begin()));
465 &*(vec_internal_stress_mean.data().begin()));
466 } else {
468 th_internal_stress, &ent, 1,
469 &*(
dataAtPts->internalStressMat->data().begin()));
471 th_actual_stress, &ent, 1,
472 &*(
dataAtPts->actualStressMat->data().begin()));
473
475 th_hydrostatic_stress, &ent, 1,
476 &*(
dataAtPts->hydrostaticStressMat->data().begin()));
478 th_deviatoric_stress, &ent, 1,
479 &*(
dataAtPts->deviatoricStressMat->data().begin()));
480 }
481
483}
#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 ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
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.
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
constexpr double t
plate stiffness
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor0IntegrationWeight()
Get integration weights.