14#include <adolc/adolc.h>
21 auto create_vec = [&]() {
22 constexpr int ghosts[] = {0};
34 return 2 * (
order - 1) + addToRule;
40 CHKERR VolumeElementForcesAndSourcesCore::preProcess();
42 if (A != PETSC_NULLPTR) {
46 if (
F != PETSC_NULLPTR) {
66 CHKERR VecAssemblyBegin(V);
74 CHKERR VolumeElementForcesAndSourcesCore::postProcess();
86 std::vector<VectorDouble> &values_at_gauss_pts,
87 std::vector<MatrixDouble> &gardient_at_gauss_pts)
90 valuesAtGaussPts(values_at_gauss_pts),
91 gradientAtGaussPts(gardient_at_gauss_pts), zeroAtType(MBVERTEX) {}
98 const int nb_base_functions = data.
getN().size2();
102 const int nb_gauss_pts = data.
getN().size1();
103 const int rank = data.
getFieldDofs()[0]->getNbOfCoeffs();
106 if (type == zeroAtType) {
107 valuesAtGaussPts.resize(nb_gauss_pts);
108 gradientAtGaussPts.resize(nb_gauss_pts);
109 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
110 valuesAtGaussPts[gg].resize(rank,
false);
111 gradientAtGaussPts[gg].resize(rank, 3,
false);
113 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
114 valuesAtGaussPts[gg].clear();
115 gradientAtGaussPts[gg].clear();
126 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
128 double &val = valuesAtGaussPts[gg][0];
130 &gradientAtGaussPts[gg](0, 1),
131 &gradientAtGaussPts[gg](0, 2));
133 for (; bb != nb_dofs; bb++) {
134 val += base_function * field_data;
135 grad(
i) += diff_base_functions(
i) * field_data;
136 ++diff_base_functions;
140 for (; bb != nb_base_functions; bb++) {
141 ++diff_base_functions;
146 }
else if (rank == 3) {
148 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
151 &valuesAtGaussPts[gg][1],
152 &valuesAtGaussPts[gg][2]);
154 &gradientAtGaussPts[gg](0, 0), &gradientAtGaussPts[gg](0, 1),
155 &gradientAtGaussPts[gg](0, 2), &gradientAtGaussPts[gg](1, 0),
156 &gradientAtGaussPts[gg](1, 1), &gradientAtGaussPts[gg](1, 2),
157 &gradientAtGaussPts[gg](2, 0), &gradientAtGaussPts[gg](2, 1),
158 &gradientAtGaussPts[gg](2, 2));
160 for (; bb != nb_dofs / 3; bb++) {
161 values(
i) += base_function * field_data(
i);
162 gradient(
i,
j) += field_data(
i) * diff_base_functions(
j);
163 ++diff_base_functions;
167 for (; bb != nb_base_functions; bb++) {
168 ++diff_base_functions;
176 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
179 for (
int dd = 0; dd < nb_dofs / rank; dd++) {
180 for (
int rr1 = 0; rr1 < rank; rr1++) {
181 valuesAtGaussPts[gg][rr1] +=
N[dd] * values[rank * dd + rr1];
182 for (
int rr2 = 0; rr2 < 3; rr2++) {
183 gradientAtGaussPts[gg](rr1, rr2) +=
184 diffN(dd, rr2) * values[rank * dd + rr1];
202 int tag,
bool jacobian,
bool ale,
206 dAta(data),
commonData(common_data),
tAg(tag), adlocReturnValue(0),
207 jAcobian(jacobian), fUnction(!jacobian), aLe(ale), fieldDisp(field_disp) {
216 CHKERR dAta.materialAdoublePtr->calculateP_PiolaKirchhoffI(
217 dAta, getNumeredEntFiniteElementPtr());
220 auto &t_P = dAta.materialAdoublePtr->t_P;
221 auto &t_invH = dAta.materialAdoublePtr->t_invH;
222 t_P(
i,
j) = t_P(
i,
k) * t_invH(
j,
k);
223 t_P(
i,
j) *= dAta.materialAdoublePtr->detH;
227 for (
int dd1 = 0; dd1 < 3; dd1++) {
228 for (
int dd2 = 0; dd2 < 3; dd2++) {
229 dAta.materialAdoublePtr->P(dd1, dd2) >>=
244 dAta.materialAdoublePtr->F.resize(3, 3,
false);
248 nbActiveVariables = 0;
249 for (
int dd1 = 0; dd1 < 3; dd1++) {
250 for (
int dd2 = 0; dd2 < 3; dd2++) {
251 dAta.materialAdoublePtr->F(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
254 dAta.materialAdoublePtr->F(dd1, dd2) += 1;
263 nbActiveVariables = 0;
265 dAta.materialAdoublePtr->h.resize(3, 3,
false);
266 for (
int dd1 = 0; dd1 < 3; dd1++) {
267 for (
int dd2 = 0; dd2 < 3; dd2++) {
268 dAta.materialAdoublePtr->h(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
273 dAta.materialAdoublePtr->H.resize(3, 3,
false);
274 for (
int dd1 = 0; dd1 < 3; dd1++) {
275 for (
int dd2 = 0; dd2 < 3; dd2++) {
276 dAta.materialAdoublePtr->H(dd1, dd2) <<= (*ptrH)[gg](dd1, dd2);
282 dAta.materialAdoublePtr->invH.resize(3, 3,
false);
284 dAta.materialAdoublePtr->detH,
285 dAta.materialAdoublePtr->invH);
287 auto &t_F = dAta.materialAdoublePtr->t_F;
288 auto &t_h = dAta.materialAdoublePtr->t_h;
289 auto &t_invH = dAta.materialAdoublePtr->t_invH;
291 t_F(
i,
j) = t_h(
i,
k) * t_invH(
k,
j);
295 CHKERR dAta.materialAdoublePtr->setUserActiveVariables(nbActiveVariables);
296 CHKERR calculateStress(gg);
312 r = ::function(
tAg, 9, nbActiveVariables, &activeVariables[0],
314 if (r < adlocReturnValue) {
316 "ADOL-C function evaluation with error r = %d", r);
322 double *jac_ptr[] = {
329 r = jacobian(
tAg, 9, nbActiveVariables, &activeVariables[0], jac_ptr);
330 if (r < adlocReturnValue) {
332 "ADOL-C function evaluation with error");
340 int row_side, EntityType row_type,
345 if (row_type != MBVERTEX)
348 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
356 dAta.materialAdoublePtr->commonDataPtr = &
commonData;
357 dAta.materialAdoublePtr->opPtr =
this;
359 int nb_gauss_pts = row_data.
getN().size1();
368 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
370 dAta.materialAdoublePtr->gG = gg;
373 if (recordTagForIntegrationPoint(gg)) {
378 if (jAcobian || (!recordTagForIntegrationPoint(gg))) {
379 activeVariables.resize(nbActiveVariables,
false);
381 for (
int dd1 = 0; dd1 < 3; dd1++) {
382 for (
int dd2 = 0; dd2 < 3; dd2++) {
383 activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
387 for (
int dd1 = 0; dd1 < 3; dd1++) {
388 for (
int dd2 = 0; dd2 < 3; dd2++) {
389 activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
392 for (
int dd1 = 0; dd1 < 3; dd1++) {
393 for (
int dd2 = 0; dd2 < 3; dd2++) {
394 activeVariables(9 + dd1 * 3 + dd2) = (*ptrH)[gg](dd1, dd2);
398 CHKERR dAta.materialAdoublePtr->setUserActiveVariables(activeVariables);
401 if (jAcobian || (!recordTagForIntegrationPoint(gg))) {
414 bool hessian,
bool ale,
bool field_disp)
417 dAta(data),
commonData(common_data),
tAg(tag), gRadient(gradient),
418 hEssian(hessian), aLe(ale), fieldDisp(field_disp) {}
423 CHKERR dAta.materialAdoublePtr->calculateElasticEnergy(
424 dAta, getNumeredEntFiniteElementPtr());
437 nbActiveVariables = 0;
438 for (
int dd1 = 0; dd1 < 3; dd1++) {
439 for (
int dd2 = 0; dd2 < 3; dd2++) {
440 dAta.materialAdoublePtr->F(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
443 dAta.materialAdoublePtr->F(dd1, dd2) += 1;
452 nbActiveVariables = 0;
454 dAta.materialAdoublePtr->h.resize(3, 3,
false);
455 for (
int dd1 = 0; dd1 < 3; dd1++) {
456 for (
int dd2 = 0; dd2 < 3; dd2++) {
457 dAta.materialAdoublePtr->h(dd1, dd2) <<= (*ptrh)[gg](dd1, dd2);
462 dAta.materialAdoublePtr->H.resize(3, 3,
false);
463 for (
int dd1 = 0; dd1 < 3; dd1++) {
464 for (
int dd2 = 0; dd2 < 3; dd2++) {
465 dAta.materialAdoublePtr->H(dd1, dd2) <<= (*ptrH)[gg](dd1, dd2);
471 dAta.materialAdoublePtr->invH.resize(3, 3,
false);
473 dAta.materialAdoublePtr->detH,
474 dAta.materialAdoublePtr->invH);
476 auto &t_F = dAta.materialAdoublePtr->t_F;
477 auto &t_h = dAta.materialAdoublePtr->t_h;
478 auto &t_invH = dAta.materialAdoublePtr->t_invH;
480 t_F(
i,
j) = t_h(
i,
k) * t_invH(
k,
j);
484 CHKERR dAta.materialAdoublePtr->setUserActiveVariables(nbActiveVariables);
498 int r = ::gradient(
tAg, nbActiveVariables, &activeVariables[0],
504 "ADOL-C function evaluation with error");
511 double *
H[nbActiveVariables];
512 for (
int n = 0;
n != nbActiveVariables;
n++) {
515 int r = ::hessian(
tAg, nbActiveVariables, &*activeVariables.begin(),
H);
520 "ADOL-C function evaluation with error");
528 int row_side, EntityType row_type,
533 if (row_type != MBVERTEX)
536 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
544 dAta.materialAdoublePtr->commonDataPtr = &
commonData;
545 dAta.materialAdoublePtr->opPtr =
this;
547 int nb_gauss_pts = row_data.
getN().size1();
556 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
558 dAta.materialAdoublePtr->gG = gg;
561 if (recordTagForIntegrationPoint(gg)) {
565 activeVariables.resize(nbActiveVariables,
false);
567 for (
int dd1 = 0; dd1 < 3; dd1++) {
568 for (
int dd2 = 0; dd2 < 3; dd2++) {
569 activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
573 for (
int dd1 = 0; dd1 < 3; dd1++) {
574 for (
int dd2 = 0; dd2 < 3; dd2++) {
575 activeVariables(dd1 * 3 + dd2) = (*ptrh)[gg](dd1, dd2);
578 for (
int dd1 = 0; dd1 < 3; dd1++) {
579 for (
int dd2 = 0; dd2 < 3; dd2++) {
580 activeVariables(9 + dd1 * 3 + dd2) = (*ptrH)[gg](dd1, dd2);
584 CHKERR dAta.materialAdoublePtr->setUserActiveVariables(activeVariables);
597 dAta(data),
commonData(common_data), aLe(false) {}
600 int row_side, EntityType row_type,
606 if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
607 iNdices.resize(nb_dofs,
false);
609 indices_ptr = &iNdices[0];
611 VectorDofs::iterator dit = dofs.begin();
612 for (
int ii = 0; dit != dofs.end(); dit++, ii++) {
613 if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
614 dAta.forcesOnlyOnEntitiesRow.end()) {
625 int row_side, EntityType row_type,
629 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
634 const int nb_dofs = row_data.
getIndices().size();
637 if ((
unsigned int)nb_dofs > 3 * row_data.
getN().size2()) {
638 SETERRQ(PETSC_COMM_SELF, 1,
"data inconsistency");
640 const int nb_base_functions = row_data.
getN().size2();
641 const int nb_gauss_pts = row_data.
getN().size1();
643 nf.resize(nb_dofs,
false);
650 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
651 double val = getVolume() * getGaussPts()(3, gg);
654 &stress(0, 0), &stress(0, 1), &stress(0, 2), &stress(1, 0),
655 &stress(1, 1), &stress(1, 2), &stress(2, 0), &stress(2, 1),
659 for (; bb != nb_dofs / 3; bb++) {
660 rhs(
i) += val * t3(
i,
j) * diff_base_functions(
j);
662 ++diff_base_functions;
664 for (; bb != nb_base_functions; bb++) {
665 ++diff_base_functions;
669 CHKERR aSemble(row_side, row_type, row_data);
681 dAta(data),
commonData(common_data), ghostVec(ghost_vec, true),
682 fieldDisp(field_disp) {}
685 int row_side, EntityType row_type,
689 if (row_type != MBVERTEX)
691 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
696 std::vector<MatrixDouble> &
F =
698 dAta.materialDoublePtr->F.resize(3, 3,
false);
702 for (
unsigned int gg = 0; gg != row_data.
getN().size1(); ++gg) {
703 double val = getVolume() * getGaussPts()(3, gg);
704 noalias(dAta.materialDoublePtr->F) =
F[gg];
706 for (
int dd = 0; dd < 3; dd++) {
707 dAta.materialDoublePtr->F(dd, dd) += 1;
710 int nb_active_variables = 0;
711 CHKERR dAta.materialDoublePtr->setUserActiveVariables(nb_active_variables);
712 CHKERR dAta.materialDoublePtr->calculateElasticEnergy(
713 dAta, getNumeredEntFiniteElementPtr());
714 energy += val * dAta.materialDoublePtr->eNergy;
717 CHKERR VecSetValue(ghostVec, 0, energy, ADD_VALUES);
726 dAta(data),
commonData(common_data), aLe(false) {}
739 const_cast<double *
>(&(col_data.
getDiffN(gg, nb_col / 3)(0, 0)));
743 &jac_stress(3 * 0 + 0, S + 0), &jac_stress(3 * 0 + 0, S + 1),
744 &jac_stress(3 * 0 + 0, S + 2), &jac_stress(3 * 0 + 1, S + 0),
745 &jac_stress(3 * 0 + 1, S + 1), &jac_stress(3 * 0 + 1, S + 2),
746 &jac_stress(3 * 0 + 2, S + 0), &jac_stress(3 * 0 + 2, S + 1),
747 &jac_stress(3 * 0 + 2, S + 2), &jac_stress(3 * 1 + 0, S + 0),
748 &jac_stress(3 * 1 + 0, S + 1), &jac_stress(3 * 1 + 0, S + 2),
749 &jac_stress(3 * 1 + 1, S + 0), &jac_stress(3 * 1 + 1, S + 1),
750 &jac_stress(3 * 1 + 1, S + 2), &jac_stress(3 * 1 + 2, S + 0),
751 &jac_stress(3 * 1 + 2, S + 1), &jac_stress(3 * 1 + 2, S + 2),
752 &jac_stress(3 * 2 + 0, S + 0), &jac_stress(3 * 2 + 0, S + 1),
753 &jac_stress(3 * 2 + 0, S + 2), &jac_stress(3 * 2 + 1, S + 0),
754 &jac_stress(3 * 2 + 1, S + 1), &jac_stress(3 * 2 + 1, S + 2),
755 &jac_stress(3 * 2 + 2, S + 0), &jac_stress(3 * 2 + 2, S + 1),
756 &jac_stress(3 * 2 + 2, S + 2));
758 &jac_stress(3 * 0 + 0, S + 3), &jac_stress(3 * 0 + 0, S + 4),
759 &jac_stress(3 * 0 + 0, S + 5), &jac_stress(3 * 0 + 1, S + 3),
760 &jac_stress(3 * 0 + 1, S + 4), &jac_stress(3 * 0 + 1, S + 5),
761 &jac_stress(3 * 0 + 2, S + 3), &jac_stress(3 * 0 + 2, S + 4),
762 &jac_stress(3 * 0 + 2, S + 5), &jac_stress(3 * 1 + 0, S + 3),
763 &jac_stress(3 * 1 + 0, S + 4), &jac_stress(3 * 1 + 0, S + 5),
764 &jac_stress(3 * 1 + 1, S + 3), &jac_stress(3 * 1 + 1, S + 4),
765 &jac_stress(3 * 1 + 1, S + 5), &jac_stress(3 * 1 + 2, S + 3),
766 &jac_stress(3 * 1 + 2, S + 4), &jac_stress(3 * 1 + 2, S + 5),
767 &jac_stress(3 * 2 + 0, S + 3), &jac_stress(3 * 2 + 0, S + 4),
768 &jac_stress(3 * 2 + 0, S + 5), &jac_stress(3 * 2 + 1, S + 3),
769 &jac_stress(3 * 2 + 1, S + 4), &jac_stress(3 * 2 + 1, S + 5),
770 &jac_stress(3 * 2 + 2, S + 3), &jac_stress(3 * 2 + 2, S + 4),
771 &jac_stress(3 * 2 + 2, S + 5));
773 &jac_stress(3 * 0 + 0, S + 6), &jac_stress(3 * 0 + 0, S + 7),
774 &jac_stress(3 * 0 + 0, S + 8), &jac_stress(3 * 0 + 1, S + 6),
775 &jac_stress(3 * 0 + 1, S + 7), &jac_stress(3 * 0 + 1, S + 8),
776 &jac_stress(3 * 0 + 2, S + 6), &jac_stress(3 * 0 + 2, S + 7),
777 &jac_stress(3 * 0 + 2, S + 8), &jac_stress(3 * 1 + 0, S + 6),
778 &jac_stress(3 * 1 + 0, S + 7), &jac_stress(3 * 1 + 0, S + 8),
779 &jac_stress(3 * 1 + 1, S + 6), &jac_stress(3 * 1 + 1, S + 7),
780 &jac_stress(3 * 1 + 1, S + 8), &jac_stress(3 * 1 + 2, S + 6),
781 &jac_stress(3 * 1 + 2, S + 7), &jac_stress(3 * 1 + 2, S + 8),
782 &jac_stress(3 * 2 + 0, S + 6), &jac_stress(3 * 2 + 0, S + 7),
783 &jac_stress(3 * 2 + 0, S + 8), &jac_stress(3 * 2 + 1, S + 6),
784 &jac_stress(3 * 2 + 1, S + 7), &jac_stress(3 * 2 + 1, S + 8),
785 &jac_stress(3 * 2 + 2, S + 6), &jac_stress(3 * 2 + 2, S + 7),
786 &jac_stress(3 * 2 + 2, S + 8));
790 &jac(0, 0), &jac(1, 0), &jac(2, 0), &jac(3, 0), &jac(4, 0), &jac(5, 0),
791 &jac(6, 0), &jac(7, 0), &jac(8, 0));
793 &jac(0, 1), &jac(1, 1), &jac(2, 1), &jac(3, 1), &jac(4, 1), &jac(5, 1),
794 &jac(6, 1), &jac(7, 1), &jac(8, 1));
796 &jac(0, 2), &jac(1, 2), &jac(2, 2), &jac(3, 2), &jac(4, 2), &jac(5, 2),
797 &jac(6, 2), &jac(7, 2), &jac(8, 2));
799 diff_ptr, &diff_ptr[1], &diff_ptr[2]);
800 for (
int dd = 0; dd != nb_col / 3; ++dd) {
801 t2_1_0(
i,
j) += t3_1_0(
i,
j,
k) * diff(
k);
802 t2_1_1(
i,
j) += t3_1_1(
i,
j,
k) * diff(
k);
803 t2_1_2(
i,
j) += t3_1_2(
i,
j,
k) * diff(
k);
818 int row_side,
int col_side, EntityType row_type, EntityType col_type,
826 int *row_indices_ptr = &row_data.
getIndices()[0];
827 int *col_indices_ptr = &col_data.
getIndices()[0];
829 if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
830 rowIndices.resize(nb_row,
false);
832 row_indices_ptr = &rowIndices[0];
834 VectorDofs::iterator dit = dofs.begin();
835 for (
int ii = 0; dit != dofs.end(); dit++, ii++) {
836 if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
837 dAta.forcesOnlyOnEntitiesRow.end()) {
843 if (!dAta.forcesOnlyOnEntitiesCol.empty()) {
844 colIndices.resize(nb_col,
false);
846 col_indices_ptr = &colIndices[0];
848 VectorDofs::iterator dit = dofs.begin();
849 for (
int ii = 0; dit != dofs.end(); dit++, ii++) {
850 if (dAta.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
851 dAta.forcesOnlyOnEntitiesCol.end()) {
858 col_indices_ptr, &
k(0, 0), ADD_VALUES);
861 if (row_side != col_side || row_type != col_type) {
866 if (!dAta.forcesOnlyOnEntitiesCol.empty()) {
867 rowIndices.resize(nb_row,
false);
869 row_indices_ptr = &rowIndices[0];
871 VectorDofs::iterator dit = dofs.begin();
872 for (
int ii = 0; dit != dofs.end(); dit++, ii++) {
873 if (dAta.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
874 dAta.forcesOnlyOnEntitiesCol.end()) {
880 if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
881 colIndices.resize(nb_col,
false);
883 col_indices_ptr = &colIndices[0];
885 VectorDofs::iterator dit = dofs.begin();
886 for (
int ii = 0; dit != dofs.end(); dit++, ii++) {
887 if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
888 dAta.forcesOnlyOnEntitiesRow.end()) {
894 trans_k.resize(nb_col, nb_row,
false);
895 noalias(trans_k) = trans(
k);
897 row_indices_ptr, &trans_k(0, 0), ADD_VALUES);
904 int row_side,
int col_side, EntityType row_type, EntityType col_type,
916 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
922 const int nb_gauss_pts = row_data.
getN().size1();
928 k.resize(nb_row, nb_col,
false);
930 jac.resize(9, nb_col,
false);
932 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
933 CHKERR getJac(col_data, gg);
934 double val = getVolume() * getGaussPts()(3, gg);
936 &jac(3 * 0 + 0, 0), &jac(3 * 0 + 0, 1), &jac(3 * 0 + 0, 2),
937 &jac(3 * 0 + 1, 0), &jac(3 * 0 + 1, 1), &jac(3 * 0 + 1, 2),
938 &jac(3 * 0 + 2, 0), &jac(3 * 0 + 2, 1), &jac(3 * 0 + 2, 2),
939 &jac(3 * 1 + 0, 0), &jac(3 * 1 + 0, 1), &jac(3 * 1 + 0, 2),
940 &jac(3 * 1 + 1, 0), &jac(3 * 1 + 1, 1), &jac(3 * 1 + 1, 2),
941 &jac(3 * 1 + 2, 0), &jac(3 * 1 + 2, 1), &jac(3 * 1 + 2, 2),
942 &jac(3 * 2 + 0, 0), &jac(3 * 2 + 0, 1), &jac(3 * 2 + 0, 2),
943 &jac(3 * 2 + 1, 0), &jac(3 * 2 + 1, 1), &jac(3 * 2 + 1, 2),
944 &jac(3 * 2 + 2, 0), &jac(3 * 2 + 2, 1), &jac(3 * 2 + 2, 2));
945 for (
int cc = 0; cc != nb_col / 3; cc++) {
948 &
k(0, 3 * cc + 0), &
k(0, 3 * cc + 1), &
k(0, 3 * cc + 2),
949 &
k(1, 3 * cc + 0), &
k(1, 3 * cc + 1), &
k(1, 3 * cc + 2),
950 &
k(2, 3 * cc + 0), &
k(2, 3 * cc + 1), &
k(2, 3 * cc + 2), 3 * nb_col);
951 for (
int rr = 0; rr != nb_row / 3; rr++) {
952 lhs(
i,
j) += val * t3_1(
i,
m,
j) * diff_base_functions(
m);
953 ++diff_base_functions;
960 CHKERR aSemble(row_side, col_side, row_type, col_type, row_data, col_data);
978 int row_side,
int col_side, EntityType row_type, EntityType col_type,
986 int *row_indices_ptr = &row_data.
getIndices()[0];
987 if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
988 rowIndices.resize(nb_row,
false);
990 row_indices_ptr = &rowIndices[0];
992 VectorDofs::iterator dit = dofs.begin();
993 for (
int ii = 0; dit != dofs.end(); dit++, ii++) {
994 if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
995 dAta.forcesOnlyOnEntitiesRow.end()) {
1001 int *col_indices_ptr = &col_data.
getIndices()[0];
1002 if (!dAta.forcesOnlyOnEntitiesCol.empty()) {
1003 colIndices.resize(nb_col,
false);
1005 col_indices_ptr = &colIndices[0];
1007 VectorDofs::iterator dit = dofs.begin();
1008 for (
int ii = 0; dit != dofs.end(); dit++, ii++) {
1009 if (dAta.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
1010 dAta.forcesOnlyOnEntitiesCol.end()) {
1011 colIndices[ii] = -1;
1025 col_indices_ptr, &
k(0, 0), ADD_VALUES);
1032 int tag,
bool jacobian,
bool ale)
1034 jacobian, ale, false) {}
1041 CHKERR dAta.materialAdoublePtr->calculatesIGma_EshelbyStress(
1042 dAta, getNumeredEntFiniteElementPtr());
1044 auto &t_sIGma = dAta.materialAdoublePtr->t_sIGma;
1045 auto &t_invH = dAta.materialAdoublePtr->t_invH;
1046 t_sIGma(
i,
j) = t_sIGma(
i,
k) * t_invH(
j,
k);
1047 t_sIGma(
i,
j) *= dAta.materialAdoublePtr->detH;
1051 for (
int dd1 = 0; dd1 < 3; dd1++) {
1052 for (
int dd2 = 0; dd2 < 3; dd2++) {
1053 dAta.materialAdoublePtr->sIGma(dd1, dd2) >>=
1089 materialAdoublePtr) {
1092 if (!materialDoublePtr) {
1094 "Pointer for materialDoublePtr not allocated");
1096 if (!materialAdoublePtr) {
1098 "Pointer for materialAdoublePtr not allocated");
1104 CHKERR it->getAttributeDataStructure(mydata);
1105 int id = it->getMeshsetId();
1112 setOfBlocks[id].materialDoublePtr = materialDoublePtr;
1113 setOfBlocks[id].materialAdoublePtr = materialAdoublePtr;
1120 const std::string element_name,
1121 const std::string spatial_position_field_name,
1122 const std::string material_position_field_name,
const bool ale) {
1127 element_name, spatial_position_field_name);
1129 element_name, spatial_position_field_name);
1131 element_name, spatial_position_field_name);
1135 element_name, material_position_field_name);
1137 element_name, material_position_field_name);
1140 element_name, material_position_field_name);
1143 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
1153 const std::string spatial_position_field_name,
1154 const std::string material_position_field_name,
const bool ale,
1155 const bool field_disp) {
1168 std::map<int, BlockData>::iterator sit =
setOfBlocks.begin();
1171 spatial_position_field_name, sit->second,
commonData,
tAg,
false, ale,
1174 spatial_position_field_name, sit->second,
commonData));
1201 spatial_position_field_name, sit->second,
commonData,
tAg,
true, ale,
1204 spatial_position_field_name, spatial_position_field_name, sit->second,
static MoFEMErrorCode calculateEnergy(DM dm, boost::shared_ptr< map< int, BlockData > > block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, SmartPetscObj< Vec > &v_energy_ptr)
ForcesAndSourcesCore::UserDataOperator UserDataOperator
static MoFEMErrorCode get_jac(EntitiesFieldData::EntData &col_data, int gg, MatrixDouble &jac_stress, MatrixDouble &jac)
Operators and data structures for non-linear elastic analysis.
#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 ...
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
@ MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
virtual bool check_field(const std::string &name) const =0
check if field is in database
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
const double n
refractive index of diffusive medium
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorShallowArrayAdaptor< double > VectorAdaptor
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Matrix adaptor.
implementation of Data Operators for Forces and Sources
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.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
constexpr auto field_name
FTensor::Index< 'm', 3 > m
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
bool sYmm
If true assume that matrix is symmetric structure.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from field data coefficients.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
const VectorInt & getIndices() const
Get global indices of dofs on entity.
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Elastic material data structure.
intrusive_ptr for managing petsc objects
FTensor::Tensor2< double *, 3, 3 > & getJac()
get element Jacobian
Volume finite element base.
data for calculation heat conductivity and heat capacity elements
common data used by volume elements
std::vector< VectorDouble > jacEnergy
std::vector< MatrixDouble > jacStress
this is simply material tangent operator
std::vector< double > eNergy
std::vector< VectorDouble > hessianEnergy
std::vector< MatrixDouble3by3 > sTress
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
Implementation of elastic (non-linear) St. Kirchhoff equation.
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode postProcess()
function is run at the end of loop
MyVolumeFE(MoFEM::Interface &m_field)
int getRule(int order)
it is used to calculate nb. of Gauss integration points
OpEnergy(const std::string field_name, BlockData &data, CommonData &common_data, SmartPetscObj< Vec > ghost_vec, bool field_disp)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpGetCommonDataAtGaussPts(const std::string field_name, CommonData &common_data)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
operator calculating deformation gradient
OpGetDataAtGaussPts(const std::string field_name, std::vector< VectorDouble > &values_at_gauss_pts, std::vector< MatrixDouble > &gradient_at_gauss_pts)
virtual MoFEMErrorCode playTag(const int gg)
Play ADOL-C tape.
OpJacobianEnergy(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool gradient, bool hessian, bool ale, bool field_disp)
virtual MoFEMErrorCode recordTag(const int gg)
Record ADOL-C tape.
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
virtual MoFEMErrorCode calculateEnergy(const int gg)
Calculate Paola-Kirchhoff I stress.
MoFEMErrorCode calculateStress(const int gg)
Calculate Paola-Kirchhoff I stress.
OpJacobianEshelbyStress(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian, bool ale)
Operator performs automatic differentiation.
virtual MoFEMErrorCode recordTag(const int gg)
Record ADOL-C tape.
virtual MoFEMErrorCode playTag(const int gg)
Play ADOL-C tape.
OpJacobianPiolaKirchhoffStress(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian, bool ale, bool field_disp)
Construct operator to calculate Piola-Kirchhoff stress or its derivatives over gradient deformation.
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
Calculate stress or jacobian at gauss points.
virtual MoFEMErrorCode calculateStress(const int gg)
Calculate Paola-Kirchhoff I stress.
OpLhsEshelby_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
OpLhsEshelby_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
OpLhsPiolaKirchhoff_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
virtual MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpLhsPiolaKirchhoff_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
OpRhsEshelbyStress(const std::string field_name, BlockData &data, CommonData &common_data)
virtual MoFEMErrorCode aSemble(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpRhsPiolaKirchhoff(const std::string field_name, BlockData &data, CommonData &common_data)
MyVolumeFE feRhs
calculate right hand side for tetrahedral elements
NonlinearElasticElement(MoFEM::Interface &m_field, short int tag)
MoFEMErrorCode addElement(const std::string element_name, const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false)
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
FTensor::Index< 'k', 3 > k
FTensor::Index< 'j', 3 > j
MoFEM::Interface & mField
MoFEMErrorCode setBlocks(boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< double > > materialDoublePtr, boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< adouble > > materialAdoublePtr)
MyVolumeFE feEnergy
calculate elastic energy
FTensor::Index< 'i', 3 > i
MoFEMErrorCode setOperators(const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false, const bool field_disp=false)
Set operators to calculate left hand tangent matrix and right hand residual.