23using namespace boost::numeric;
41 nF.resize(nb_dofs,
false);
45 const int nb_gauss_pts = data.
getN().size1();
54 const double normal_stiffness =
commonDataPtr->springStiffnessNormal;
55 const double tangent_stiffness =
commonDataPtr->springStiffnessTangent;
64 auto t_solution_at_gauss_point =
66 auto t_init_solution_at_gauss_point =
70 auto t_normal_1 = getFTensor1FromMat<3>(
commonDataPtr->normalVector);
72 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
73 t_normal(
i) = t_normal_1(
i);
75 const double normal_length = std::sqrt(t_normal(
k) * t_normal(
k));
76 t_normal(
i) = t_normal(
i) / normal_length;
77 t_normal_projection(
i,
j) = t_normal(
i) * t_normal(
j);
78 t_tangent_projection(
i,
j) = t_normal_projection(
i,
j);
79 t_tangent_projection(0, 0) -= 1;
80 t_tangent_projection(1, 1) -= 1;
81 t_tangent_projection(2, 2) -= 1;
84 t_displacement_at_gauss_point(
i) =
85 t_solution_at_gauss_point(
i) - t_init_solution_at_gauss_point(
i);
87 t_displacement_at_gauss_point(
i) = t_solution_at_gauss_point(
i);
90 double w = t_w * normal_length;
102 for (
int rr = 0; rr != nb_dofs / 3; ++rr) {
104 t_nf(
i) +=
w * t_base_func * normal_stiffness *
105 t_normal_projection(
i,
j) * t_displacement_at_gauss_point(
j);
106 t_nf(
i) -=
w * t_base_func * tangent_stiffness *
107 t_tangent_projection(
i,
j) * t_displacement_at_gauss_point(
j);
117 ++t_solution_at_gauss_point;
118 ++t_init_solution_at_gauss_point;
135 const int row_nb_dofs = row_data.
getIndices().size();
139 const int col_nb_dofs = col_data.
getIndices().size();
143 if (dAta.tRis.find(getFEEntityHandle()) == dAta.tRis.end()) {
147 CHKERR commonDataPtr->getBlockData(dAta);
149 locKs.resize(row_nb_dofs, col_nb_dofs,
false);
153 const int row_nb_gauss_pts = row_data.
getN().size1();
154 if (!row_nb_gauss_pts)
158 auto t_w = getFTensor0IntegrationWeight();
169 auto t_tangent1_ptr = getFTensor1Tangent1AtGaussPts();
171 t_tangent1(
i) = t_tangent1_ptr(
i);
177 auto t_tangent2_ptr = getFTensor1Tangent2AtGaussPts();
181 auto normal_original_slave =
182 getFTensor1FromMat<3>(commonDataPtr->normalVector);
187 for (
int gg = 0; gg != row_nb_gauss_pts; gg++) {
188 t_normal(
i) = normal_original_slave(
i);
189 const double normal_length = std::sqrt(t_normal(
k) * t_normal(
k));
190 t_normal(
i) = t_normal(
i) / normal_length;
191 t_normal_projection(
i,
j) = t_normal(
i) * t_normal(
j);
192 t_tangent_projection(
i,
j) = t_normal_projection(
i,
j);
193 t_tangent_projection(0, 0) -= 1;
194 t_tangent_projection(1, 1) -= 1;
195 t_tangent_projection(2, 2) -= 1;
198 double w = t_w * normal_length;
199 if (getNumeredEntFiniteElementPtr()->getEntType() == MBTRI) {
205 for (
int rr = 0; rr != row_nb_dofs / 3; rr++) {
206 auto assemble_m = getFTensor2FromArray<3, 3, 3>(locKs, 3 * rr);
208 for (
int cc = 0; cc != col_nb_dofs / 3; cc++) {
209 assemble_m(
i,
j) +=
w * t_row_base_func * t_col_base_func *
210 commonDataPtr->springStiffnessNormal *
211 t_normal_projection(
i,
j);
212 assemble_m(
i,
j) -=
w * t_row_base_func * t_col_base_func *
213 commonDataPtr->springStiffnessTangent *
214 t_tangent_projection(
i,
j);
222 ++normal_original_slave;
229 if (row_side != col_side || row_type != col_type) {
230 transLocKs.resize(col_nb_dofs, row_nb_dofs,
false);
231 noalias(transLocKs) = trans(locKs);
248 const int row_nb_dofs = row_data.
getIndices().size();
252 const int col_nb_dofs = col_data.
getIndices().size();
256 if (dAta.tRis.find(getFEEntityHandle()) == dAta.tRis.end()) {
260 CHKERR commonDataPtr->getBlockData(dAta);
262 locKs.resize(row_nb_dofs, col_nb_dofs,
false);
266 const int row_nb_gauss_pts = row_data.
getN().size1();
267 if (!row_nb_gauss_pts)
271 auto t_w = getFTensor0IntegrationWeight();
278 auto make_vec_der = [&](
auto t_N,
auto t_1,
auto t_2) {
286 auto make_vec_der_2 = [&](
auto t_N,
auto t_1,
auto t_2) {
289 const double normal_norm = std::sqrt(t_normal(
i) * t_normal(
i));
297 t_n(
i,
j) / normal_norm - t_normal(
i) * t_n(
k,
j) * t_normal(
k) /
298 (normal_norm * normal_norm * normal_norm);
302 const double normal_stiffness = commonDataPtr->springStiffnessNormal;
303 const double tangent_stiffness = commonDataPtr->springStiffnessTangent;
313 auto t_solution_at_gauss_point =
314 getFTensor1FromMat<3>(*commonDataPtr->xAtPts);
315 auto t_init_solution_at_gauss_point =
316 getFTensor1FromMat<3>(*commonDataPtr->xInitAtPts);
318 auto t_1 = getFTensor1FromMat<3>(commonDataPtr->tangent1);
319 auto t_2 = getFTensor1FromMat<3>(commonDataPtr->tangent2);
323 auto normal_at_gp = getFTensor1FromMat<3>(commonDataPtr->normalVector);
326 for (
int gg = 0; gg != row_nb_gauss_pts; gg++) {
328 t_normal(
i) = normal_at_gp(
i);
329 const double normal_length = std::sqrt(t_normal(
k) * t_normal(
k));
330 t_normal(
i) = t_normal(
i) / normal_length;
332 t_normal_projection(
i,
j) = t_normal(
i) * t_normal(
j);
333 t_tangent_projection(
i,
j) = t_normal_projection(
i,
j);
334 t_tangent_projection(0, 0) -= 1;
335 t_tangent_projection(1, 1) -= 1;
336 t_tangent_projection(2, 2) -= 1;
338 double w = 0.5 * t_w;
340 t_displacement_at_gauss_point(
i) =
341 t_solution_at_gauss_point(
i) - t_init_solution_at_gauss_point(
i);
345 for (
int rr = 0; rr != row_nb_dofs / 3; rr++) {
348 auto assemble_m = getFTensor2FromArray<3, 3, 3>(locKs, 3 * rr);
349 for (
int cc = 0; cc != col_nb_dofs / 3; cc++) {
350 auto t_d_n = make_vec_der(t_N, t_1, t_2);
351 auto t_d_n_2 = make_vec_der_2(t_N, t_1, t_2);
353 assemble_m(
i,
j) -=
w * normal_length * t_col_base_func *
354 normal_stiffness * t_row_base_func *
355 t_normal_projection(
i,
j);
357 assemble_m(
i,
j) +=
w * normal_length * t_col_base_func *
358 tangent_stiffness * t_row_base_func *
359 t_tangent_projection(
i,
j);
362 w * (normal_stiffness - tangent_stiffness) * t_row_base_func *
363 (t_d_n(
i,
j) * (t_normal(
k) * t_displacement_at_gauss_point(
k)) +
364 normal_length * t_normal(
i) *
365 (t_d_n_2(
k,
j) * t_displacement_at_gauss_point(
k)));
367 assemble_m(
i,
j) +=
w * (t_d_n(
l,
j) * t_normal(
l)) *
368 tangent_stiffness * t_row_base_func *
t_kd(
i,
k) *
369 t_displacement_at_gauss_point(
k);
379 ++t_solution_at_gauss_point;
380 ++t_init_solution_at_gauss_point;
398 if (dataAtSpringPts->faceRowData ==
nullptr)
401 if (row_type != MBVERTEX)
404 row_nb_dofs = dataAtSpringPts->faceRowData->getIndices().size();
411 nb_gauss_pts = dataAtSpringPts->faceRowData->getN().size1();
413 nb_base_fun_row = dataAtSpringPts->faceRowData->getFieldData().size() / 3;
416 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col,
false);
420 CHKERR iNtegrate(*(dataAtSpringPts->faceRowData), col_data);
423 CHKERR aSsemble(*(dataAtSpringPts->faceRowData), col_data);
435 const int *row_indices = &*row_data.
getIndices().data().begin();
437 const int *col_indices = &*col_data.
getIndices().data().begin();
439 auto &data = *dataAtSpringPts;
440 if (!data.forcesOnlyOnEntitiesRow.empty()) {
441 rowIndices.resize(row_nb_dofs,
false);
443 row_indices = &rowIndices[0];
445 VectorDofs::iterator dit = dofs.begin();
446 for (
int ii = 0; dit != dofs.end(); ++dit, ++ii) {
447 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
448 data.forcesOnlyOnEntitiesRow.end()) {
454 if (!data.forcesOnlyOnEntitiesCol.empty()) {
455 colIndices.resize(col_nb_dofs,
false);
457 col_indices = &colIndices[0];
459 VectorDofs::iterator dit = dofs.begin();
460 for (
int ii = 0; dit != dofs.end(); ++dit, ++ii) {
461 if (data.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
462 data.forcesOnlyOnEntitiesCol.end()) {
468 Mat
B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
469 : getFEMethod()->snes_B;
472 &*NN.data().begin(), ADD_VALUES);
489 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2), &
m(
r + 1,
c + 0),
490 &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2), &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1),
497 const double normal_stiffness = dataAtSpringPts->springStiffnessNormal;
498 const double tangent_stiffness = dataAtSpringPts->springStiffnessTangent;
501 auto t_solution_at_gauss_point =
502 getFTensor1FromMat<3>(*dataAtSpringPts->xAtPts);
503 auto t_init_solution_at_gauss_point =
504 getFTensor1FromMat<3>(*dataAtSpringPts->xInitAtPts);
507 auto t_normal_1 = getFTensor1FromMat<3>(dataAtSpringPts->normalVector);
508 auto t_w = getFTensor0IntegrationWeight();
510 auto t_inv_H = getFTensor2FromMat<3, 3>(*dataAtSpringPts->invHMat);
511 auto t_F = getFTensor2FromMat<3, 3>(*dataAtSpringPts->FMat);
512 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
514 t_normal(
i) = t_normal_1(
i);
516 const double normal_length = std::sqrt(t_normal(
k) * t_normal(
k));
517 t_normal(
i) = t_normal(
i) / normal_length;
518 t_normal_projection(
i,
j) = t_normal(
i) * t_normal(
j);
519 t_tangent_projection(
i,
j) = t_normal_projection(
i,
j);
520 t_tangent_projection(0, 0) -= 1;
521 t_tangent_projection(1, 1) -= 1;
522 t_tangent_projection(2, 2) -= 1;
524 t_displacement_at_gauss_point(
i) =
525 t_solution_at_gauss_point(
i) - t_init_solution_at_gauss_point(
i);
527 double w = 0.5 * t_w * normal_length;
532 for (; bbc != nb_base_fun_col; ++bbc) {
537 for (; bbr != nb_base_fun_row; ++bbr) {
539 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
541 t_assemble(
i,
j) -=
w * t_row_base * t_col_base * normal_stiffness *
542 t_F(
k,
i) * t_normal_projection(
k,
j);
543 t_assemble(
i,
j) +=
w * t_row_base * t_col_base * tangent_stiffness *
544 t_F(
k,
i) * t_tangent_projection(
k,
j);
546 t_assemble(
i,
j) -=
w * normal_stiffness * t_row_base * t_inv_H(
k,
i) *
547 t_col_diff_base(
k) * t_normal_projection(
j,
l) *
548 t_displacement_at_gauss_point(
l);
549 t_assemble(
i,
j) +=
w * tangent_stiffness * t_row_base * t_inv_H(
k,
i) *
550 t_col_diff_base(
k) * t_tangent_projection(
j,
l) *
551 t_displacement_at_gauss_point(
l);
559 ++t_solution_at_gauss_point;
560 ++t_init_solution_at_gauss_point;
576 const int *row_indices = &*row_data.
getIndices().data().begin();
578 const int *col_indices = &*col_data.
getIndices().data().begin();
580 auto &data = *dataAtSpringPts;
581 if (!data.forcesOnlyOnEntitiesRow.empty()) {
582 rowIndices.resize(row_nb_dofs,
false);
584 row_indices = &rowIndices[0];
586 VectorDofs::iterator dit = dofs.begin();
587 for (
int ii = 0; dit != dofs.end(); ++dit, ++ii) {
588 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
589 data.forcesOnlyOnEntitiesRow.end()) {
595 if (!data.forcesOnlyOnEntitiesCol.empty()) {
596 colIndices.resize(col_nb_dofs,
false);
598 col_indices = &colIndices[0];
600 VectorDofs::iterator dit = dofs.begin();
601 for (
int ii = 0; dit != dofs.end(); ++dit, ++ii) {
602 if (data.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
603 data.forcesOnlyOnEntitiesCol.end()) {
609 Mat
B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
610 : getFEMethod()->snes_B;
613 &*NN.data().begin(), ADD_VALUES);
630 nb_gauss_pts = row_data.
getN().size1();
635 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col,
false);
638 diagonal_block = (row_type == col_type) && (row_side == col_side);
640 if (col_type == MBVERTEX) {
641 dataAtSpringPts->faceRowData = &row_data;
642 CHKERR loopSideVolumes(sideFeName, *sideFe);
646 CHKERR iNtegrate(row_data, col_data);
649 CHKERR aSsemble(row_data, col_data);
667 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2), &
m(
r + 1,
c + 0),
668 &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2), &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1),
672 auto make_vec_der = [&](
auto t_N,
auto t_1,
auto t_2) {
680 auto make_vec_der_2 = [&](
auto t_N,
auto t_1,
auto t_2) {
683 const double normal_norm = std::sqrt(t_normal(
i) * t_normal(
i));
691 t_n(
i,
j) / normal_norm - t_normal(
i) * t_n(
k,
j) * t_normal(
k) /
692 (normal_norm * normal_norm * normal_norm);
696 dataAtSpringPts->faceRowData =
nullptr;
697 CHKERR loopSideVolumes(sideFeName, *sideFe);
699 auto t_F = getFTensor2FromMat<3, 3>(*dataAtSpringPts->FMat);
701 auto t_w = getFTensor0IntegrationWeight();
702 auto t_1 = getFTensor1FromMat<3>(dataAtSpringPts->tangent1);
703 auto t_2 = getFTensor1FromMat<3>(dataAtSpringPts->tangent2);
705 const double normal_stiffness = dataAtSpringPts->springStiffnessNormal;
706 const double tangent_stiffness = dataAtSpringPts->springStiffnessTangent;
715 auto t_solution_at_gauss_point =
716 getFTensor1FromMat<3>(*dataAtSpringPts->xAtPts);
717 auto t_init_solution_at_gauss_point =
718 getFTensor1FromMat<3>(*dataAtSpringPts->xInitAtPts);
721 auto t_normal_1 = getFTensor1FromMat<3>(dataAtSpringPts->normalVector);
724 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
725 t_normal(
i) = t_normal_1(
i);
727 const double normal_length = std::sqrt(t_normal(
k) * t_normal(
k));
728 t_normal(
i) = t_normal(
i) / normal_length;
729 t_normal_projection(
i,
j) = t_normal(
i) * t_normal(
j);
730 t_tangent_projection(
i,
j) = t_normal_projection(
i,
j);
731 t_tangent_projection(0, 0) -= 1;
732 t_tangent_projection(1, 1) -= 1;
733 t_tangent_projection(2, 2) -= 1;
735 t_displacement_at_gauss_point(
i) =
736 t_solution_at_gauss_point(
i) - t_init_solution_at_gauss_point(
i);
737 double val = 0.5 * t_w;
742 for (; bbc != nb_base_fun_col; ++bbc) {
747 for (; bbr != nb_base_fun_row; ++bbr) {
749 auto t_d_n = make_vec_der(t_N, t_1, t_2);
751 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
753 auto t_d_n_2 = make_vec_der_2(t_N, t_1, t_2);
755 t_assemble(
i,
j) += val * normal_length * t_col_base_func *
756 normal_stiffness * t_row_base_func * t_F(
k,
i) *
757 t_normal_projection(
k,
j);
759 t_assemble(
i,
j) -= val * normal_length * t_col_base_func *
760 tangent_stiffness * t_row_base_func * t_F(
k,
i) *
761 t_tangent_projection(
k,
j);
764 val * (normal_stiffness - tangent_stiffness) * t_row_base_func *
766 (t_d_n(
l,
j) * (t_normal(
k) * t_displacement_at_gauss_point(
k)) +
767 normal_length * t_normal(
l) *
768 (t_d_n_2(
k,
j) * t_displacement_at_gauss_point(
k)));
770 t_assemble(
i,
j) -= val * (t_d_n(
l,
j) * t_normal(
l)) *
771 tangent_stiffness * t_row_base_func * t_F(
k,
i) *
772 t_displacement_at_gauss_point(
k);
784 ++t_solution_at_gauss_point;
785 ++t_init_solution_at_gauss_point;
805 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2), &
m(
r + 1,
c + 0),
806 &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2), &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1),
810 const double normal_stiffness = dataAtSpringPts->springStiffnessNormal;
811 const double tangent_stiffness = dataAtSpringPts->springStiffnessTangent;
820 auto t_solution_at_gauss_point =
821 getFTensor1FromMat<3>(*dataAtSpringPts->xAtPts);
822 auto t_init_solution_at_gauss_point =
823 getFTensor1FromMat<3>(*dataAtSpringPts->xInitAtPts);
826 auto t_normal_1 = getFTensor1FromMat<3>(dataAtSpringPts->normalVector);
828 auto t_w = getFTensor0IntegrationWeight();
830 auto t_h = getFTensor2FromMat<3, 3>(*dataAtSpringPts->hMat);
831 auto t_inv_H = getFTensor2FromMat<3, 3>(*dataAtSpringPts->invHMat);
835 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
836 t_normal(
i) = t_normal_1(
i);
838 const double normal_length = std::sqrt(t_normal(
k) * t_normal(
k));
839 t_normal(
i) = t_normal(
i) / normal_length;
840 t_normal_projection(
i,
j) = t_normal(
i) * t_normal(
j);
841 t_tangent_projection(
i,
j) = t_normal_projection(
i,
j);
842 t_tangent_projection(0, 0) -= 1;
843 t_tangent_projection(1, 1) -= 1;
844 t_tangent_projection(2, 2) -= 1;
846 t_displacement_at_gauss_point(
i) =
847 t_solution_at_gauss_point(
i) - t_init_solution_at_gauss_point(
i);
849 double val = 0.5 * t_w * normal_length;
854 for (; bbc != nb_base_fun_col; ++bbc) {
859 for (; bbr != nb_base_fun_row; ++bbr) {
861 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
864 val * t_row_base * normal_stiffness * t_inv_H(
l,
j) *
865 t_col_diff_base(
m) * t_inv_H(
m,
i) * t_h(
k,
l) *
866 (t_normal_projection(
k, q) * t_displacement_at_gauss_point(q));
869 val * t_row_base * tangent_stiffness * t_inv_H(
l,
j) *
870 t_col_diff_base(
m) * t_inv_H(
m,
i) * t_h(
k,
l) *
871 (t_tangent_projection(
k, q) * t_displacement_at_gauss_point(q));
880 ++t_solution_at_gauss_point;
881 ++t_init_solution_at_gauss_point;
898 const int nb_integration_pts = getGaussPts().size2();
900 auto t_h = getFTensor2FromMat<3, 3>(*commonDataPtr->hMat);
901 auto t_H = getFTensor2FromMat<3, 3>(*commonDataPtr->HMat);
903 commonDataPtr->detHVec->resize(nb_integration_pts,
false);
904 commonDataPtr->invHMat->resize(9, nb_integration_pts,
false);
905 commonDataPtr->FMat->resize(9, nb_integration_pts,
false);
907 commonDataPtr->detHVec->clear();
908 commonDataPtr->invHMat->clear();
909 commonDataPtr->FMat->clear();
912 auto t_invH = getFTensor2FromMat<3, 3>(*commonDataPtr->invHMat);
913 auto t_F = getFTensor2FromMat<3, 3>(*commonDataPtr->FMat);
915 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
918 t_F(
i,
j) = t_h(
i,
k) * t_invH(
k,
j);
935 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
946 nF.resize(nbRows,
false);
950 nbIntegrationPts = getGaussPts().size2();
953 CHKERR iNtegrate(row_data);
956 CHKERR aSsemble(row_data);
964 CHKERR loopSideVolumes(sideFeName, *sideFe);
971 CHKERR dataAtPts->getBlockData(dAta);
974 auto t_w = getFTensor0IntegrationWeight();
981 const double normal_stiffness = dataAtPts->springStiffnessNormal;
982 const double tangent_stiffness = dataAtPts->springStiffnessTangent;
991 auto t_solution_at_gauss_point = getFTensor1FromMat<3>(*dataAtPts->xAtPts);
992 auto t_init_solution_at_gauss_point =
993 getFTensor1FromMat<3>(*dataAtPts->xInitAtPts);
996 auto t_normal_1 = getFTensor1FromMat<3>(dataAtPts->normalVector);
998 auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
1001 for (
int gg = 0; gg != nbIntegrationPts; ++gg) {
1002 t_normal(
i) = t_normal_1(
i);
1004 const double normal_length = std::sqrt(t_normal(
k) * t_normal(
k));
1005 t_normal(
i) = t_normal(
i) / normal_length;
1006 t_normal_projection(
i,
j) = t_normal(
i) * t_normal(
j);
1007 t_tangent_projection(
i,
j) = t_normal_projection(
i,
j);
1008 t_tangent_projection(0, 0) -= 1;
1009 t_tangent_projection(1, 1) -= 1;
1010 t_tangent_projection(2, 2) -= 1;
1012 t_displacement_at_gauss_point(
i) =
1013 t_solution_at_gauss_point(
i) - t_init_solution_at_gauss_point(
i);
1015 double w = t_w * 0.5 * normal_length;
1024 for (
int rr = 0; rr != nb_dofs / 3; ++rr) {
1026 t_nf(
i) -=
w * t_base_func * normal_stiffness * t_F(
k,
i) *
1027 t_normal_projection(
k,
j) * t_displacement_at_gauss_point(
j);
1028 t_nf(
i) +=
w * t_base_func * tangent_stiffness * t_F(
k,
i) *
1029 t_tangent_projection(
k,
j) * t_displacement_at_gauss_point(
j);
1039 ++t_solution_at_gauss_point;
1040 ++t_init_solution_at_gauss_point;
1055 const int *row_indices = &*row_data.
getIndices().data().begin();
1057 auto &data = *dataAtPts;
1058 if (!data.forcesOnlyOnEntitiesRow.empty()) {
1059 rowIndices.resize(nbRows,
false);
1061 row_indices = &rowIndices[0];
1063 VectorDofs::iterator dit = dofs.begin();
1064 for (
int ii = 0; dit != dofs.end(); ++dit, ++ii) {
1065 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
1066 data.forcesOnlyOnEntitiesRow.end()) {
1067 rowIndices[ii] = -1;
1072 CHKERR VecSetOption(getSNESf(), VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
1085 PetscFunctionReturn(0);
1087 ngp = data.
getN().size1();
1091 if (
type == MBVERTEX) {
1092 dataAtIntegrationPts->tangent1.resize(3, ngp,
false);
1093 dataAtIntegrationPts->tangent1.clear();
1095 dataAtIntegrationPts->tangent2.resize(3, ngp,
false);
1096 dataAtIntegrationPts->tangent2.clear();
1099 auto t_1 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent1);
1100 auto t_2 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent2);
1102 for (
unsigned int gg = 0; gg != ngp; ++gg) {
1106 for (
unsigned int dd = 0;
dd != nb_dofs; ++
dd) {
1107 t_1(
i) += t_dof(
i) * t_N(0);
1108 t_2(
i) += t_dof(
i) * t_N(1);
1126 ngp = data.
getN().size1();
1132 if (
type == MBVERTEX) {
1133 dataAtIntegrationPts->normalVector.resize(3, ngp,
false);
1134 dataAtIntegrationPts->normalVector.clear();
1137 auto normal_original_slave =
1138 getFTensor1FromMat<3>(dataAtIntegrationPts->normalVector);
1140 auto tangent_1 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent1);
1141 auto tangent_2 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent2);
1143 for (
unsigned int gg = 0; gg != ngp; ++gg) {
1144 normal_original_slave(
i) =
1146 ++normal_original_slave;
1160 if (col_type != MBVERTEX)
1163 dataAtSpringPts->faceRowData = &row_data;
1164 CHKERR loopSideVolumes(sideFeName, *sideFe);
1172 const std::string mesh_nodals_positions) {
1183 mesh_nodals_positions);
1188 if (
bit->getName().compare(0, 9,
"SPRING_BC") == 0) {
1199 const std::string mesh_nodals_positions, Range &spring_triangles) {
1209 mesh_nodals_positions);
1211 mesh_nodals_positions);
1213 mesh_nodals_positions);
1223 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr,
1224 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr,
1225 const std::string
field_name,
const std::string mesh_nodals_positions,
double stiffness_scale) {
1230 boost::shared_ptr<MetaSpringBC::DataAtIntegrationPtsSprings> commonDataPtr =
1231 boost::make_shared<MetaSpringBC::DataAtIntegrationPtsSprings>(m_field, stiffness_scale);
1232 CHKERR commonDataPtr->getParameters();
1241 for (
auto &sitSpring : commonDataPtr->mapSpring) {
1243 fe_spring_lhs_ptr->getOpPtrVector().push_back(
1245 fe_spring_lhs_ptr->getOpPtrVector().push_back(
1247 fe_spring_lhs_ptr->getOpPtrVector().push_back(
1250 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1252 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1254 commonDataPtr->xInitAtPts));
1255 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1257 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1259 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1268 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr_dx,
1269 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr_dX,
1270 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr,
1271 boost::shared_ptr<MetaSpringBC::DataAtIntegrationPtsSprings>
1272 data_at_integration_pts,
1273 const std::string
field_name,
const std::string mesh_nodals_positions,
1274 std::string side_fe_name) {
1279 CHKERR data_at_integration_pts->getParameters();
1290 for (
auto &sitSpring : data_at_integration_pts->mapSpring) {
1292 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideRhs =
1293 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(m_field);
1295 feMatSideRhs->getOpPtrVector().push_back(
1297 mesh_nodals_positions, data_at_integration_pts->HMat));
1298 feMatSideRhs->getOpPtrVector().push_back(
1303 mesh_nodals_positions, data_at_integration_pts));
1305 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1307 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1310 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1312 data_at_integration_pts->xAtPts));
1313 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1315 mesh_nodals_positions, data_at_integration_pts->xInitAtPts));
1317 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1319 feMatSideRhs, side_fe_name, sitSpring.second));
1321 fe_spring_lhs_ptr_dx->getOpPtrVector().push_back(
1323 data_at_integration_pts->xAtPts));
1324 fe_spring_lhs_ptr_dx->getOpPtrVector().push_back(
1326 mesh_nodals_positions, data_at_integration_pts->xInitAtPts));
1327 fe_spring_lhs_ptr_dx->getOpPtrVector().push_back(
1329 fe_spring_lhs_ptr_dx->getOpPtrVector().push_back(
1332 fe_spring_lhs_ptr_dx->getOpPtrVector().push_back(
1334 mesh_nodals_positions));
1336 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideLhs_dx =
1337 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(m_field);
1339 feMatSideLhs_dx->getOpPtrVector().push_back(
1341 mesh_nodals_positions, data_at_integration_pts->HMat));
1342 feMatSideLhs_dx->getOpPtrVector().push_back(
1347 mesh_nodals_positions, data_at_integration_pts));
1349 feMatSideLhs_dx->getOpPtrVector().push_back(
1351 mesh_nodals_positions,
field_name, data_at_integration_pts));
1353 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1355 data_at_integration_pts->xAtPts));
1356 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1358 mesh_nodals_positions, data_at_integration_pts->xInitAtPts));
1359 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1361 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1364 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1366 data_at_integration_pts,
1367 feMatSideLhs_dx, side_fe_name));
1369 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideLhs_dX =
1370 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(m_field);
1372 feMatSideLhs_dX->getOpPtrVector().push_back(
1374 mesh_nodals_positions, data_at_integration_pts->HMat));
1375 feMatSideLhs_dX->getOpPtrVector().push_back(
1380 mesh_nodals_positions, data_at_integration_pts));
1382 feMatSideLhs_dX->getOpPtrVector().push_back(
1384 mesh_nodals_positions,
1385 data_at_integration_pts));
1387 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1389 mesh_nodals_positions, mesh_nodals_positions,
1390 data_at_integration_pts, feMatSideLhs_dX, side_fe_name));
Kronecker Delta class symmetric.
#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.
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 add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
add entities to finite element
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_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
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 bool check_field(const std::string &name) const =0
check if field is in database
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const FTensor::Tensor2< T, Dim, Dim > Vec
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
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)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
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 getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
constexpr auto VecSetValues
constexpr auto MatSetValues
const double r
rate factor
constexpr auto field_name
FTensor::Index< 'm', 3 > m
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 & 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 filed data coeffinects.
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
const VectorInt & getIndices() const
Get global indices of dofs on entity.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor0IntegrationWeight()
Get integration weights.
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.