19 using namespace MoFEM;
22 HookeElement::OpCalculateStrainAle::OpCalculateStrainAle(
23 const std::string row_field,
const std::string col_field,
24 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
26 dataAtPts(data_at_pts) {
27 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE],
false);
30 MoFEMErrorCode HookeElement::OpCalculateStrainAle::doWork(
int row_side,
38 const int nb_integration_pts = getGaussPts().size2();
39 auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
40 auto t_H = getFTensor2FromMat<3, 3>(*dataAtPts->HMat);
42 dataAtPts->detHVec->resize(nb_integration_pts,
false);
43 dataAtPts->invHMat->resize(9, nb_integration_pts,
false);
44 dataAtPts->FMat->resize(9, nb_integration_pts,
false);
45 dataAtPts->smallStrainMat->resize(6, nb_integration_pts,
false);
48 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
49 auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
50 auto t_strain = getFTensor2SymmetricFromMat<3>(*dataAtPts->smallStrainMat);
52 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
55 t_F(
i,
j) = t_h(
i,
k) * t_invH(
k,
j);
56 t_strain(
i,
j) = (t_F(
i,
j) || t_F(
j,
i)) / 2.;
72 HookeElement::OpCalculateEnergy::OpCalculateEnergy(
73 const std::string row_field,
const std::string col_field,
74 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
77 dataAtPts(data_at_pts), ghostVec(ghost_vec, true) {
78 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE],
false);
81 MoFEMErrorCode HookeElement::OpCalculateEnergy::doWork(
int row_side,
87 const int nb_integration_pts = getGaussPts().size2();
88 auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
89 auto t_cauchy_stress =
90 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
91 dataAtPts->energyVec->resize(nb_integration_pts,
false);
93 &*(dataAtPts->energyVec->data().begin()));
98 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
99 t_energy = (t_strain(
i,
j) * t_cauchy_stress(
i,
j)) / 2.;
105 if (ghostVec.get()) {
107 double vol = getVolume();
109 auto t_w = getFTensor0IntegrationWeight();
110 auto &det_H = *dataAtPts->detHVec;
112 &*(dataAtPts->energyVec->data().begin()));
114 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
116 double a = t_w * vol;
120 energy +=
a * t_energy;
124 CHKERR VecSetValue(ghostVec, 0, energy, ADD_VALUES);
130 HookeElement::OpCalculateEshelbyStress::OpCalculateEshelbyStress(
131 const std::string row_field,
const std::string col_field,
132 boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
134 dataAtPts(data_at_pts) {
135 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE],
false);
139 int row_side, EntityType row_type,
EntData &row_data) {
142 const int nb_integration_pts = getGaussPts().size2();
144 &*(dataAtPts->energyVec->data().begin()));
145 auto t_cauchy_stress =
146 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
147 auto t_F = getFTensor2FromMat<3, 3>(*(dataAtPts->FMat));
148 dataAtPts->eshelbyStressMat->resize(9, nb_integration_pts,
false);
149 auto t_eshelby_stress =
150 getFTensor2FromMat<3, 3>(*(dataAtPts->eshelbyStressMat));
156 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
157 t_eshelby_stress(
i,
j) = -t_F(
k,
i) * t_cauchy_stress(
k,
j);
158 t_eshelby_stress(0, 0) += t_energy;
159 t_eshelby_stress(1, 1) += t_energy;
160 t_eshelby_stress(2, 2) += t_energy;
169 HookeElement::OpAssemble::OpAssemble(
170 const std::string row_field,
const std::string col_field,
171 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
const char type,
174 dataAtPts(data_at_pts) {}
176 MoFEMErrorCode HookeElement::OpAssemble::doWork(
int row_side,
int col_side,
200 K.resize(nbRows, nbCols,
false);
204 nbIntegrationPts = getGaussPts().size2();
206 if (row_side == col_side && row_type == col_type) {
213 CHKERR iNtegrate(row_data, col_data);
216 CHKERR aSsemble(row_data, col_data);
232 nF.resize(nbRows,
false);
236 nbIntegrationPts = getGaussPts().size2();
239 CHKERR iNtegrate(row_data);
242 CHKERR aSsemble(row_data);
265 const int *row_indices = &*row_data.
getIndices().data().begin();
267 const int *col_indices = &*col_data.
getIndices().data().begin();
269 auto &data = *dataAtPts;
270 if (!data.forcesOnlyOnEntitiesRow.empty()) {
271 rowIndices.resize(nbRows,
false);
273 row_indices = &rowIndices[0];
275 VectorDofs::iterator dit = dofs.begin();
276 for (
int ii = 0; dit != dofs.end(); dit++, ii++) {
277 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
278 data.forcesOnlyOnEntitiesRow.end()) {
284 if (!data.forcesOnlyOnEntitiesCol.empty()) {
285 colIndices.resize(nbCols,
false);
287 col_indices = &colIndices[0];
289 VectorDofs::iterator dit = dofs.begin();
290 for (
int ii = 0; dit != dofs.end(); dit++, ii++) {
291 if (data.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
292 data.forcesOnlyOnEntitiesCol.end()) {
298 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
299 : getFEMethod()->snes_B;
302 &*
K.data().begin(), ADD_VALUES);
304 if (!isDiag && sYmm) {
307 transK.resize(
K.size2(),
K.size1(),
false);
308 noalias(transK) = trans(
K);
310 &*transK.data().begin(), ADD_VALUES);
319 const int *row_indices = &*row_data.
getIndices().data().begin();
321 auto &data = *dataAtPts;
322 if (!data.forcesOnlyOnEntitiesRow.empty()) {
323 rowIndices.resize(nbRows,
false);
325 row_indices = &rowIndices[0];
327 VectorDofs::iterator dit = dofs.begin();
328 for (
int ii = 0; dit != dofs.end(); dit++, ii++) {
329 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
330 data.forcesOnlyOnEntitiesRow.end()) {
336 Vec F = getFEMethod()->snes_f;
342 HookeElement::OpRhs_dx::OpRhs_dx(
343 const std::string row_field,
const std::string col_field,
344 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
345 :
OpAssemble(row_field, col_field, data_at_pts, OPROW) {}
352 &
v(
r + 0), &
v(
r + 1), &
v(
r + 2));
361 double vol = getVolume();
363 auto t_w = getFTensor0IntegrationWeight();
367 const int row_nb_base_fun = row_data.
getN().size2();
368 auto t_cauchy_stress =
369 getFTensor2SymmetricFromMat<3>(*dataAtPts->cauchyStressMat);
372 for (
int gg = 0; gg != nbIntegrationPts; ++gg) {
375 double a = t_w * vol;
376 auto t_nf = get_tensor1(nF, 0);
379 for (; rr != nbRows / 3; ++rr) {
380 t_nf(
i) +=
a * t_row_diff_base(
j) * t_cauchy_stress(
i,
j);
385 for (; rr != row_nb_base_fun; ++rr)
395 HookeElement::OpAleRhs_dx::OpAleRhs_dx(
396 const std::string row_field,
const std::string col_field,
397 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
398 :
OpAssemble(row_field, col_field, data_at_pts, OPROW) {}
405 &
v(
r + 0), &
v(
r + 1), &
v(
r + 2));
412 double vol = getVolume();
414 auto t_w = getFTensor0IntegrationWeight();
418 const int row_nb_base_fun = row_data.
getN().size2();
419 auto t_cauchy_stress =
420 getFTensor2SymmetricFromMat<3>(*dataAtPts->cauchyStressMat);
421 auto &det_H = *dataAtPts->detHVec;
422 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
425 for (
int gg = 0; gg != nbIntegrationPts; ++gg) {
428 double a = t_w * vol * det_H[gg];
429 auto t_nf = get_tensor1(nF, 0);
432 for (; rr != nbRows / 3; ++rr) {
434 t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
435 t_nf(
i) +=
a * t_row_diff_base_pulled(
j) * t_cauchy_stress(
i,
j);
440 for (; rr != row_nb_base_fun; ++rr)
451 HookeElement::OpAleRhs_dX::OpAleRhs_dX(
452 const std::string row_field,
const std::string col_field,
453 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
454 :
OpAssemble(row_field, col_field, data_at_pts, OPROW) {}
461 &
v(
r + 0), &
v(
r + 1), &
v(
r + 2));
468 double vol = getVolume();
470 auto t_w = getFTensor0IntegrationWeight();
474 const int row_nb_base_fun = row_data.
getN().size2();
475 auto t_eshelby_stress =
476 getFTensor2FromMat<3, 3>(*dataAtPts->eshelbyStressMat);
477 auto &det_H = *dataAtPts->detHVec;
478 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
481 for (
int gg = 0; gg != nbIntegrationPts; ++gg) {
484 double a = t_w * vol * det_H[gg];
485 auto t_nf = get_tensor1(nF, 0);
488 for (; rr != nbRows / 3; ++rr) {
490 t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
491 t_nf(
i) +=
a * t_row_diff_base_pulled(
j) * t_eshelby_stress(
i,
j);
496 for (; rr != row_nb_base_fun; ++rr)
509 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr) {
514 "Pointer to block of sets is null");
519 CHKERR it->getAttributeDataStructure(mydata);
520 int id = it->getMeshsetId();
521 auto &block_data = (*block_sets_ptr)[id];
524 block_data.tEts,
true);
526 block_data.E = mydata.
data.Young;
527 block_data.PoissonRatio = mydata.
data.Poisson;
535 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
536 const std::string element_name,
const std::string x_field,
537 const std::string X_field,
const bool ale) {
542 "Pointer to block of sets is null");
556 for (
auto &
m : (*block_sets_ptr)) {
565 boost::shared_ptr<ForcesAndSourcesCore> fe_lhs_ptr,
566 boost::shared_ptr<ForcesAndSourcesCore> fe_rhs_ptr,
567 boost::shared_ptr<map<int, BlockData>> block_sets_ptr,
568 const std::string x_field,
const std::string X_field,
const bool ale,
569 const bool field_disp,
const EntityType
type,
570 boost::shared_ptr<DataAtIntegrationPts> data_at_pts) {
575 "Pointer to block of sets is null");
578 data_at_pts = boost::make_shared<DataAtIntegrationPts>();
581 if (ale == PETSC_FALSE) {
582 if (
type == MBPRISM) {
583 boost::shared_ptr<MatrixDouble> inv_jac_ptr(
new MatrixDouble);
584 fe_lhs_ptr->getOpPtrVector().push_back(
586 fe_lhs_ptr->getOpPtrVector().push_back(
589 fe_lhs_ptr->getOpPtrVector().push_back(
591 x_field, x_field, block_sets_ptr, data_at_pts));
592 fe_lhs_ptr->getOpPtrVector().push_back(
595 if (
type == MBPRISM) {
596 boost::shared_ptr<MatrixDouble> inv_jac_ptr(
new MatrixDouble);
597 fe_lhs_ptr->getOpPtrVector().push_back(
599 fe_lhs_ptr->getOpPtrVector().push_back(
602 fe_lhs_ptr->getOpPtrVector().push_back(
604 fe_lhs_ptr->getOpPtrVector().push_back(
606 x_field, x_field, block_sets_ptr, data_at_pts));
607 fe_lhs_ptr->getOpPtrVector().push_back(
609 fe_lhs_ptr->getOpPtrVector().push_back(
611 fe_lhs_ptr->getOpPtrVector().push_back(
613 fe_lhs_ptr->getOpPtrVector().push_back(
615 fe_lhs_ptr->getOpPtrVector().push_back(
617 fe_lhs_ptr->getOpPtrVector().push_back(
619 fe_lhs_ptr->getOpPtrVector().push_back(
621 fe_lhs_ptr->getOpPtrVector().push_back(
623 fe_lhs_ptr->getOpPtrVector().push_back(
625 fe_lhs_ptr->getOpPtrVector().push_back(
632 if (ale == PETSC_FALSE) {
633 if (
type == MBPRISM) {
634 boost::shared_ptr<MatrixDouble> inv_jac_ptr(
new MatrixDouble);
635 fe_rhs_ptr->getOpPtrVector().push_back(
637 fe_rhs_ptr->getOpPtrVector().push_back(
640 fe_rhs_ptr->getOpPtrVector().push_back(
642 fe_rhs_ptr->getOpPtrVector().push_back(
644 block_sets_ptr, data_at_pts));
646 fe_rhs_ptr->getOpPtrVector().push_back(
649 fe_rhs_ptr->getOpPtrVector().push_back(
652 fe_rhs_ptr->getOpPtrVector().push_back(
654 fe_rhs_ptr->getOpPtrVector().push_back(
655 new OpRhs_dx(x_field, x_field, data_at_pts));
657 if (
type == MBPRISM) {
658 boost::shared_ptr<MatrixDouble> inv_jac_ptr(
new MatrixDouble);
659 fe_rhs_ptr->getOpPtrVector().push_back(
661 fe_rhs_ptr->getOpPtrVector().push_back(
664 fe_rhs_ptr->getOpPtrVector().push_back(
666 fe_rhs_ptr->getOpPtrVector().push_back(
668 block_sets_ptr, data_at_pts));
669 fe_rhs_ptr->getOpPtrVector().push_back(
671 fe_rhs_ptr->getOpPtrVector().push_back(
673 fe_rhs_ptr->getOpPtrVector().push_back(
675 fe_rhs_ptr->getOpPtrVector().push_back(
677 fe_rhs_ptr->getOpPtrVector().push_back(
679 fe_rhs_ptr->getOpPtrVector().push_back(
681 fe_rhs_ptr->getOpPtrVector().push_back(
690 DM dm, boost::shared_ptr<map<int, BlockData>> block_sets_ptr,
691 const std::string x_field,
const std::string X_field,
const bool ale,
700 boost::shared_ptr<DataAtIntegrationPts> data_at_pts(
701 new DataAtIntegrationPts());
704 boost::make_shared<VolumeElementForcesAndSourcesCore>(*m_field_ptr);
705 fe_ptr->getRuleHook = [](
const double,
const double,
const double o) {
710 "MESH_NODE_POSITIONS", *fe_ptr,
true,
false,
false,
false);
714 using FatPrismElementForcesAndSourcesCore::
715 FatPrismElementForcesAndSourcesCore;
716 int getRuleTrianglesOnly(
int order) {
return 2 *
order; }
717 int getRuleThroughThickness(
int order) {
return 2 *
order; }
720 boost::shared_ptr<ForcesAndSourcesCore> prism_fe_ptr(
723 auto push_ops = [&](boost::shared_ptr<ForcesAndSourcesCore> fe_ptr,
726 boost::shared_ptr<MatrixDouble> inv_jac_ptr(
new MatrixDouble);
727 if (ale == PETSC_FALSE) {
728 if (
type == MBPRISM) {
729 fe_ptr->getOpPtrVector().push_back(
731 fe_ptr->getOpPtrVector().push_back(
734 fe_ptr->getOpPtrVector().push_back(
737 x_field, x_field, block_sets_ptr, data_at_pts));
739 fe_ptr->getOpPtrVector().push_back(
742 fe_ptr->getOpPtrVector().push_back(
745 fe_ptr->getOpPtrVector().push_back(
747 fe_ptr->getOpPtrVector().push_back(
750 if (
type == MBPRISM) {
751 fe_ptr->getOpPtrVector().push_back(
753 fe_ptr->getOpPtrVector().push_back(
756 fe_ptr->getOpPtrVector().push_back(
759 x_field, x_field, block_sets_ptr, data_at_pts));
760 fe_ptr->getOpPtrVector().push_back(
762 fe_ptr->getOpPtrVector().push_back(
764 fe_ptr->getOpPtrVector().push_back(
766 fe_ptr->getOpPtrVector().push_back(
772 CHKERR push_ops(fe_ptr, MBTET);
773 CHKERR push_ops(prism_fe_ptr, MBPRISM);
775 CHKERR VecZeroEntries(v_energy);
777 fe_ptr->snes_ctx = SnesMethod::CTX_SNESNONE;
781 CHKERR VecAssemblyBegin(v_energy);
782 CHKERR VecAssemblyEnd(v_energy);
795 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2), &
m(
r + 1,
c + 0),
796 &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2), &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1),
808 double vol = getVolume();
811 auto t_w = getFTensor0IntegrationWeight();
815 const int row_nb_base_fun = row_data.
getN().size2();
817 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
818 auto &det_H = *dataAtPts->detHVec;
820 auto get_eshelby_stress_dx = [
this]() {
824 for (
int ii = 0; ii != 3; ++ii)
825 for (
int jj = 0; jj != 3; ++jj)
826 for (
int kk = 0; kk != 3; ++kk)
827 for (
int ll = 0; ll != 3; ++ll)
828 t_eshelby_stress_dx.ptr(ii, jj, kk, ll) =
829 &(*dataAtPts->eshelbyStress_dx)(mm++, 0);
830 return t_eshelby_stress_dx;
833 auto t_eshelby_stress_dx = get_eshelby_stress_dx();
836 for (
int gg = 0; gg != nbIntegrationPts; ++gg) {
839 double a = t_w * vol * det_H[gg];
843 for (; rr != nbRows / 3; ++rr) {
846 auto t_m = get_tensor2(
K, 3 * rr, 0);
849 t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
852 t_row_stress_dx(
i,
k,
l) =
853 a * t_row_diff_base_pulled(
j) * t_eshelby_stress_dx(
i,
j,
k,
l);
859 for (
int cc = 0; cc != nbCols / 3; ++cc) {
861 t_m(
i,
k) += t_row_stress_dx(
i,
k,
l) * t_col_diff_base(
l);
874 for (; rr != row_nb_base_fun; ++rr)
879 ++t_eshelby_stress_dx;
885 HookeElement::OpAleLhsWithDensity_dX_dX::OpAleLhsWithDensity_dX_dX(
886 const std::string row_field,
const std::string col_field,
887 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
888 boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
889 boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts,
const double rho_n,
891 :
OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false),
892 rhoAtGaussPtsPtr(rho_at_gauss_pts),
893 rhoGradAtGaussPtsPtr(rho_grad_at_gauss_pts), rhoN(rho_n), rHo0(rho_0) {}
896 HookeElement::OpAleLhsWithDensity_dX_dX::iNtegrate(
EntData &row_data,
904 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2), &
m(
r + 1,
c + 0),
905 &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2), &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1),
915 double vol = getVolume();
918 auto t_w = getFTensor0IntegrationWeight();
922 const int row_nb_base_fun = row_data.
getN().size2();
927 auto t_grad_rho = getFTensor1FromMat<3>(*rhoGradAtGaussPtsPtr);
929 auto t_eshelby_stress =
930 getFTensor2FromMat<3, 3>(*dataAtPts->eshelbyStressMat);
932 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
933 auto &det_H = *dataAtPts->detHVec;
936 for (
int gg = 0; gg != nbIntegrationPts; ++gg) {
939 double a = t_w * vol * det_H[gg];
941 const double stress_dho_coef = (rhoN /
rho);
947 for (; rr != nbRows / 3; ++rr) {
950 auto t_m = get_tensor2(
K, 3 * rr, 0);
953 t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
956 t_row_stress(
i) =
a * t_row_diff_base_pulled(
j) * t_eshelby_stress(
i,
j);
963 for (
int cc = 0; cc != nbCols / 3; ++cc) {
966 t_row_stress(
i) * stress_dho_coef * t_grad_rho(
k) * t_col_base;
978 for (; rr != row_nb_base_fun; ++rr)
992 HookeElement::OpAleLhsWithDensity_dx_dX::OpAleLhsWithDensity_dx_dX(
993 const std::string row_field,
const std::string col_field,
994 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
995 boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
996 boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts,
const double rho_n,
998 :
OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false),
999 rhoAtGaussPtsPtr(rho_at_gauss_pts),
1000 rhoGradAtGaussPtsPtr(rho_grad_at_gauss_pts), rhoN(rho_n), rHo0(rho_0) {}
1003 HookeElement::OpAleLhsWithDensity_dx_dX::iNtegrate(
EntData &row_data,
1011 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2), &
m(
r + 1,
c + 0),
1012 &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2), &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1),
1022 double vol = getVolume();
1025 auto t_w = getFTensor0IntegrationWeight();
1029 const int row_nb_base_fun = row_data.
getN().size2();
1032 auto t_grad_rho = getFTensor1FromMat<3>(*rhoGradAtGaussPtsPtr);
1033 auto t_cauchy_stress =
1034 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1036 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1037 auto &det_H = *dataAtPts->detHVec;
1040 for (
int gg = 0; gg != nbIntegrationPts; ++gg) {
1043 double a = t_w * vol * det_H[gg];
1045 const double stress_dho_coef = (rhoN /
rho);
1049 for (; rr != nbRows / 3; ++rr) {
1052 auto t_m = get_tensor2(
K, 3 * rr, 0);
1055 t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
1058 t_row_stress(
i) =
a * t_row_diff_base_pulled(
j) * t_cauchy_stress(
i,
j);
1063 for (
int cc = 0; cc != nbCols / 3; ++cc) {
1066 t_row_stress(
i) * stress_dho_coef * t_grad_rho(
k) * t_col_base;
1078 for (; rr != row_nb_base_fun; ++rr)
1094 HookeElement::OpCalculateStiffnessScaledByDensityField::
1095 OpCalculateStiffnessScaledByDensityField(
1096 const std::string row_field,
const std::string col_field,
1097 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
1098 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1099 boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
const double rho_n,
1103 blockSetsPtr(block_sets_ptr), dataAtPts(data_at_pts),
1104 rhoAtGaussPtsPtr(rho_at_gauss_pts), rhoN(rho_n), rHo0(rho_0) {
1105 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE],
false);
1108 MoFEMErrorCode HookeElement::OpCalculateStiffnessScaledByDensityField::doWork(
1109 int row_side, EntityType row_type,
EntData &row_data) {
1112 if (!rhoAtGaussPtsPtr)
1113 SETERRQ(PETSC_COMM_SELF, 1,
"Calculate density with MWLS first.");
1115 for (
auto &
m : (*blockSetsPtr)) {
1117 if (
m.second.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1118 m.second.tEts.end()) {
1122 const int nb_integration_pts = getGaussPts().size2();
1123 dataAtPts->stiffnessMat->resize(36, nb_integration_pts,
false);
1127 const double young =
m.second.E;
1128 const double poisson =
m.second.PoissonRatio;
1133 const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
1140 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1142 t_D(
i,
j,
k,
l) = 0.;
1144 t_D(0, 0, 0, 0) = 1 - poisson;
1145 t_D(1, 1, 1, 1) = 1 - poisson;
1146 t_D(2, 2, 2, 2) = 1 - poisson;
1148 t_D(0, 1, 0, 1) = 0.5 * (1 - 2 * poisson);
1149 t_D(0, 2, 0, 2) = 0.5 * (1 - 2 * poisson);
1150 t_D(1, 2, 1, 2) = 0.5 * (1 - 2 * poisson);
1152 t_D(0, 0, 1, 1) = poisson;
1153 t_D(1, 1, 0, 0) = poisson;
1154 t_D(0, 0, 2, 2) = poisson;
1155 t_D(2, 2, 0, 0) = poisson;
1156 t_D(1, 1, 2, 2) = poisson;
1157 t_D(2, 2, 1, 1) = poisson;
1160 t_D(
i,
j,
k,
l) *= coefficient * pow(
rho / rHo0, rhoN);