24 #ifndef __HOOKE_INTERNAL_STRESS_ELEMENT_HPP__
25 #define __HOOKE_INTERNAL_STRESS_ELEMENT_HPP__
57 const std::string col_field,
58 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
77 boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
90 typedef boost::function<
102 const std::string row_field,
const std::string col_field,
103 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
128 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
129 map<int, BlockData> &block_sets_ptr,
131 bool save_mean =
false,
bool is_ale =
false,
132 bool is_field_disp =
true)
141 template <
class ELEMENT>
152 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
153 map<int, BlockData> &block_sets_ptr,
155 std::vector<EntityHandle> &map_gauss_pts,
156 bool is_ale =
false,
bool is_field_disp =
true)
168 int row_side, EntityType row_type,
EntData &row_data) {
171 const int nb_integration_pts = row_data.
getN().size1();
175 const int val_num = 9 * nb_integration_pts;
176 std::vector<double> def_vals(val_num, 0.0);
177 Tag th_internal_stress;
180 MB_TAG_CREAT | MB_TAG_SPARSE, &*def_vals.begin());
182 dataAtPts->internalStressMat->resize(9, nb_integration_pts,
false);
185 th_internal_stress, &ent, 1,
186 &*(
dataAtPts->internalStressMat->data().begin()));
201 &
v(
r + 0), &
v(
r + 1), &
v(
r + 2));
204 const int nb_integration_pts = getGaussPts().size2();
206 auto t_internal_stress =
207 getFTensor2FromMat<3, 3>(*dataAtPts->internalStressMat);
210 double vol = getVolume();
211 auto t_w = getFTensor0IntegrationWeight();
213 nF.resize(nbRows,
false);
218 const int row_nb_base_fun = row_data.
getN().size2();
220 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
223 double a = t_w * vol;
224 auto t_nf = get_tensor1(nF, 0);
227 for (; rr != nbRows / 3; ++rr) {
228 t_nf(
i) +=
a * t_row_diff_base(
j) * t_internal_stress(
i,
j);
233 for (; rr != row_nb_base_fun; ++rr)
246 int row_side, EntityType row_type,
EntData &row_data) {
253 auto tensor_to_tensor = [](
const auto &t1,
auto &t2) {
257 t2(0, 1) = t2(1, 0) = t1(1, 0);
258 t2(0, 2) = t2(2, 0) = t1(2, 0);
259 t2(1, 2) = t2(2, 1) = t1(2, 1);
262 const int nb_integration_pts = getGaussPts().size2();
263 auto t_coords = getFTensor1CoordsAtGaussPts();
270 dataAtPts->internalStressMat->resize(9, nb_integration_pts,
false);
271 auto t_internal_stress =
272 getFTensor2FromMat<3, 3>(*dataAtPts->internalStressMat);
276 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
278 auto t_fun_strain = strainFun(t_coords);
279 t_internal_stress_symm(
i,
j) = -t_D(
i,
j,
k,
l) * t_fun_strain(
k,
l);
281 tensor_to_tensor(t_internal_stress_symm, t_internal_stress);
292 int row_side, EntityType row_type,
EntData &row_data) {
295 if (row_type != MBVERTEX) {
301 auto tensor_to_tensor = [](
const auto &t1,
auto &t2) {
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);
310 auto tensor_to_vector = [](
const auto &
t,
auto &
v) {
319 auto get_tag_handle = [&](
auto name,
auto size) {
321 std::vector<double> def_vals(size, 0.0);
322 CHKERR outputMesh.tag_get_handle(name, size, MB_TYPE_DOUBLE,
th,
323 MB_TAG_CREAT | MB_TAG_SPARSE,
328 const int nb_integration_pts = row_data.
getN().size1();
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);
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);
347 auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
348 auto t_H = getFTensor2FromMat<3, 3>(*dataAtPts->HMat);
350 dataAtPts->stiffnessMat->resize(36, 1,
false);
353 for (
auto &
m : (blockSetsPtr)) {
354 const double young =
m.second.E;
355 const double poisson =
m.second.PoissonRatio;
357 const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
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;
379 auto t_internal_stress =
380 getFTensor2FromMat<3, 3>(*dataAtPts->internalStressMat);
382 dataAtPts->actualStressMat->resize(9, nb_integration_pts,
false);
383 auto t_actual_stress = getFTensor2FromMat<3, 3>(*dataAtPts->actualStressMat);
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);
395 t_internal_stress_mean(
i,
j) = 0.;
396 t_actual_stress_mean(
i,
j) = 0.;
398 auto t_w = getFTensor0IntegrationWeight();
400 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
403 t_small_strain_symm(
i,
j) = (t_h(
i,
j) || t_h(
j,
i)) / 2.;
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.;
412 if (isALE || !isFieldDisp) {
413 for (
auto ii : {0, 1, 2}) {
414 t_small_strain_symm(ii, ii) -= 1.;
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);
423 t_actual_stress(
i,
j) *= scaleFactor;
424 t_internal_stress(
i,
j) *= scaleFactor;
426 double hydrostatic_pressure =
427 (t_actual_stress(0, 0) + t_actual_stress(1, 1) +
428 t_actual_stress(2, 2)) /
431 t_hydrostatic_stress(
i,
j) = 0.;
432 for (
auto ii : {0, 1, 2}) {
433 t_hydrostatic_stress(ii, ii) += hydrostatic_pressure;
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;
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);
452 vec_actual_stress_mean.resize(9,
false);
453 vec_actual_stress_mean.clear();
456 vec_internal_stress_mean.resize(9,
false);
457 vec_internal_stress_mean.clear();
459 tensor_to_vector(t_actual_stress_mean, vec_actual_stress_mean);
460 tensor_to_vector(t_internal_stress_mean, vec_internal_stress_mean);
462 CHKERR outputMesh.tag_set_data(th_actual_stress_mean, &ent, 1,
463 &*(vec_actual_stress_mean.data().begin()));
464 CHKERR outputMesh.tag_set_data(th_internal_stress_mean, &ent, 1,
465 &*(vec_internal_stress_mean.data().begin()));
467 CHKERR outputMesh.tag_set_data(
468 th_internal_stress, &ent, 1,
469 &*(dataAtPts->internalStressMat->data().begin()));
470 CHKERR outputMesh.tag_set_data(
471 th_actual_stress, &ent, 1,
472 &*(dataAtPts->actualStressMat->data().begin()));
474 CHKERR outputMesh.tag_set_data(
475 th_hydrostatic_stress, &ent, 1,
476 &*(dataAtPts->hydrostaticStressMat->data().begin()));
477 CHKERR outputMesh.tag_set_data(
478 th_deviatoric_stress, &ent, 1,
479 &*(dataAtPts->deviatoricStressMat->data().begin()));
485 template <
class ELEMENT>
491 if (
type != MBVERTEX) {
495 auto tensor_to_tensor = [](
const auto &t1,
auto &t2) {
499 t2(0, 1) = t2(1, 0) = t1(1, 0);
500 t2(0, 2) = t2(2, 0) = t1(2, 0);
501 t2(1, 2) = t2(2, 1) = t1(2, 1);
504 auto get_tag_handle = [&](
auto name,
auto size) {
506 std::vector<double> def_vals(size, 0.0);
507 CHKERR postProcMesh.tag_get_handle(name, size, MB_TYPE_DOUBLE,
th,
508 MB_TAG_CREAT | MB_TAG_SPARSE,
513 auto th_stress = get_tag_handle(
"STRESS", 9);
514 auto th_strain = get_tag_handle(
"STRAIN", 9);
515 auto th_psi = get_tag_handle(
"ENERGY", 1);
516 auto th_disp = get_tag_handle(
"DISPLACEMENT", 3);
518 const int nb_integration_pts = mapGaussPts.size();
525 auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
526 auto t_H = getFTensor2FromMat<3, 3>(*dataAtPts->HMat);
528 dataAtPts->stiffnessMat->resize(36, 1,
false);
531 for (
auto &
m : (blockSetsPtr)) {
532 const double young =
m.second.E;
533 const double poisson =
m.second.PoissonRatio;
535 const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
537 t_D(
i,
j,
k,
l) = 0.;
538 t_D(0, 0, 0, 0) = t_D(1, 1, 1, 1) = t_D(2, 2, 2, 2) = 1 - poisson;
539 t_D(0, 1, 0, 1) = t_D(0, 2, 0, 2) = t_D(1, 2, 1, 2) =
540 0.5 * (1 - 2 * poisson);
541 t_D(0, 0, 1, 1) = t_D(1, 1, 0, 0) = t_D(0, 0, 2, 2) = t_D(2, 2, 0, 0) =
542 t_D(1, 1, 2, 2) = t_D(2, 2, 1, 1) = poisson;
543 t_D(
i,
j,
k,
l) *= coefficient;
557 auto t_spat_pos = getFTensor1FromMat<3>(*dataAtPts->spatPosMat);
558 auto t_mesh_node_pos = getFTensor1FromMat<3>(*dataAtPts->meshNodePosMat);
561 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
570 t_small_strain_symm(
i,
j) = (t_h(
i,
j) || t_h(
j,
i)) / 2.;
574 t_F(
i,
j) = t_h(
i,
k) * t_invH(
k,
j);
575 t_small_strain_symm(
i,
j) = (t_F(
i,
j) || t_F(
j,
i)) / 2.;
579 t_small_strain_symm(0, 0) -= 1;
580 t_small_strain_symm(1, 1) -= 1;
581 t_small_strain_symm(2, 2) -= 1;
584 t_stress_symm(
i,
j) = t_D(
i,
j,
k,
l) * t_small_strain_symm(
k,
l);
585 tensor_to_tensor(t_stress_symm, t_stress);
586 tensor_to_tensor(t_small_strain_symm, t_small_strain);
588 t_disp(
i) = t_spat_pos(
i) - t_mesh_node_pos(
i);
590 const double psi = 0.5 * t_stress_symm(
i,
j) * t_small_strain_symm(
i,
j);
592 CHKERR postProcMesh.tag_set_data(th_psi, &mapGaussPts[gg], 1, &psi);
593 CHKERR postProcMesh.tag_set_data(th_stress, &mapGaussPts[gg], 1,
595 CHKERR postProcMesh.tag_set_data(th_strain, &mapGaussPts[gg], 1,
596 &t_small_strain(0, 0));
597 CHKERR postProcMesh.tag_set_data(th_disp, &mapGaussPts[gg], 1, &t_disp(0));