11using namespace boost::numeric;
29 nF.resize(nb_dofs,
false);
33 const int nb_gauss_pts = data.
getN().size1();
42 const double normal_stiffness =
commonDataPtr->springStiffnessNormal;
43 const double tangent_stiffness =
commonDataPtr->springStiffnessTangent;
52 auto t_solution_at_gauss_point =
54 auto t_init_solution_at_gauss_point =
58 auto t_normal_1 = getFTensor1FromMat<3>(
commonDataPtr->normalVector);
60 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
61 t_normal(
i) = t_normal_1(
i);
63 const double normal_length = std::sqrt(t_normal(
k) * t_normal(
k));
64 t_normal(
i) = t_normal(
i) / normal_length;
65 t_normal_projection(
i,
j) = t_normal(
i) * t_normal(
j);
66 t_tangent_projection(
i,
j) = t_normal_projection(
i,
j);
67 t_tangent_projection(0, 0) -= 1;
68 t_tangent_projection(1, 1) -= 1;
69 t_tangent_projection(2, 2) -= 1;
72 t_displacement_at_gauss_point(
i) =
73 t_solution_at_gauss_point(
i) - t_init_solution_at_gauss_point(
i);
75 t_displacement_at_gauss_point(
i) = t_solution_at_gauss_point(
i);
78 double w = t_w * normal_length;
90 for (
int rr = 0; rr != nb_dofs / 3; ++rr) {
92 t_nf(
i) += w * t_base_func * normal_stiffness *
93 t_normal_projection(
i,
j) * t_displacement_at_gauss_point(
j);
94 t_nf(
i) -= w * t_base_func * tangent_stiffness *
95 t_tangent_projection(
i,
j) * t_displacement_at_gauss_point(
j);
105 ++t_solution_at_gauss_point;
106 ++t_init_solution_at_gauss_point;
123 const int row_nb_dofs = row_data.
getIndices().size();
127 const int col_nb_dofs = col_data.
getIndices().size();
131 if (dAta.tRis.find(getFEEntityHandle()) == dAta.tRis.end()) {
135 CHKERR commonDataPtr->getBlockData(dAta);
137 locKs.resize(row_nb_dofs, col_nb_dofs,
false);
141 const int row_nb_gauss_pts = row_data.
getN().size1();
142 if (!row_nb_gauss_pts)
146 auto t_w = getFTensor0IntegrationWeight();
157 auto t_tangent1_ptr = getFTensor1Tangent1AtGaussPts();
159 t_tangent1(
i) = t_tangent1_ptr(
i);
165 auto t_tangent2_ptr = getFTensor1Tangent2AtGaussPts();
169 auto normal_original_slave =
170 getFTensor1FromMat<3>(commonDataPtr->normalVector);
175 for (
int gg = 0; gg != row_nb_gauss_pts; gg++) {
176 t_normal(
i) = normal_original_slave(
i);
177 const double normal_length = std::sqrt(t_normal(
k) * t_normal(
k));
178 t_normal(
i) = t_normal(
i) / normal_length;
179 t_normal_projection(
i,
j) = t_normal(
i) * t_normal(
j);
180 t_tangent_projection(
i,
j) = t_normal_projection(
i,
j);
181 t_tangent_projection(0, 0) -= 1;
182 t_tangent_projection(1, 1) -= 1;
183 t_tangent_projection(2, 2) -= 1;
186 double w = t_w * normal_length;
187 if (getNumeredEntFiniteElementPtr()->getEntType() == MBTRI) {
193 for (
int rr = 0; rr != row_nb_dofs / 3; rr++) {
194 auto assemble_m = getFTensor2FromArray<3, 3, 3>(locKs, 3 * rr);
196 for (
int cc = 0; cc != col_nb_dofs / 3; cc++) {
197 assemble_m(
i,
j) += w * t_row_base_func * t_col_base_func *
198 commonDataPtr->springStiffnessNormal *
199 t_normal_projection(
i,
j);
200 assemble_m(
i,
j) -= w * t_row_base_func * t_col_base_func *
201 commonDataPtr->springStiffnessTangent *
202 t_tangent_projection(
i,
j);
210 ++normal_original_slave;
217 if (row_side != col_side || row_type != col_type) {
218 transLocKs.resize(col_nb_dofs, row_nb_dofs,
false);
219 noalias(transLocKs) = trans(locKs);
236 const int row_nb_dofs = row_data.
getIndices().size();
240 const int col_nb_dofs = col_data.
getIndices().size();
244 if (dAta.tRis.find(getFEEntityHandle()) == dAta.tRis.end()) {
248 CHKERR commonDataPtr->getBlockData(dAta);
250 locKs.resize(row_nb_dofs, col_nb_dofs,
false);
254 const int row_nb_gauss_pts = row_data.
getN().size1();
255 if (!row_nb_gauss_pts)
259 auto t_w = getFTensor0IntegrationWeight();
266 auto make_vec_der = [&](
auto t_N,
auto t_1,
auto t_2) {
274 auto make_vec_der_2 = [&](
auto t_N,
auto t_1,
auto t_2) {
277 const double normal_norm = std::sqrt(t_normal(
i) * t_normal(
i));
285 t_n(
i,
j) / normal_norm - t_normal(
i) * t_n(
k,
j) * t_normal(
k) /
286 (normal_norm * normal_norm * normal_norm);
290 const double normal_stiffness = commonDataPtr->springStiffnessNormal;
291 const double tangent_stiffness = commonDataPtr->springStiffnessTangent;
301 auto t_solution_at_gauss_point =
302 getFTensor1FromMat<3>(*commonDataPtr->xAtPts);
303 auto t_init_solution_at_gauss_point =
304 getFTensor1FromMat<3>(*commonDataPtr->xInitAtPts);
306 auto t_1 = getFTensor1FromMat<3>(commonDataPtr->tangent1);
307 auto t_2 = getFTensor1FromMat<3>(commonDataPtr->tangent2);
311 auto normal_at_gp = getFTensor1FromMat<3>(commonDataPtr->normalVector);
314 for (
int gg = 0; gg != row_nb_gauss_pts; gg++) {
316 t_normal(
i) = normal_at_gp(
i);
317 const double normal_length = std::sqrt(t_normal(
k) * t_normal(
k));
318 t_normal(
i) = t_normal(
i) / normal_length;
320 t_normal_projection(
i,
j) = t_normal(
i) * t_normal(
j);
321 t_tangent_projection(
i,
j) = t_normal_projection(
i,
j);
322 t_tangent_projection(0, 0) -= 1;
323 t_tangent_projection(1, 1) -= 1;
324 t_tangent_projection(2, 2) -= 1;
326 double w = 0.5 * t_w;
328 t_displacement_at_gauss_point(
i) =
329 t_solution_at_gauss_point(
i) - t_init_solution_at_gauss_point(
i);
333 for (
int rr = 0; rr != row_nb_dofs / 3; rr++) {
336 auto assemble_m = getFTensor2FromArray<3, 3, 3>(locKs, 3 * rr);
337 for (
int cc = 0; cc != col_nb_dofs / 3; cc++) {
338 auto t_d_n = make_vec_der(t_N, t_1, t_2);
339 auto t_d_n_2 = make_vec_der_2(t_N, t_1, t_2);
341 assemble_m(
i,
j) -= w * normal_length * t_col_base_func *
342 normal_stiffness * t_row_base_func *
343 t_normal_projection(
i,
j);
345 assemble_m(
i,
j) += w * normal_length * t_col_base_func *
346 tangent_stiffness * t_row_base_func *
347 t_tangent_projection(
i,
j);
350 w * (normal_stiffness - tangent_stiffness) * t_row_base_func *
351 (t_d_n(
i,
j) * (t_normal(
k) * t_displacement_at_gauss_point(
k)) +
352 normal_length * t_normal(
i) *
353 (t_d_n_2(
k,
j) * t_displacement_at_gauss_point(
k)));
355 assemble_m(
i,
j) += w * (t_d_n(
l,
j) * t_normal(
l)) *
356 tangent_stiffness * t_row_base_func *
t_kd(
i,
k) *
357 t_displacement_at_gauss_point(
k);
367 ++t_solution_at_gauss_point;
368 ++t_init_solution_at_gauss_point;
386 if (dataAtSpringPts->faceRowData ==
nullptr)
389 if (row_type != MBVERTEX)
392 row_nb_dofs = dataAtSpringPts->faceRowData->getIndices().size();
399 nb_gauss_pts = dataAtSpringPts->faceRowData->getN().size1();
401 nb_base_fun_row = dataAtSpringPts->faceRowData->getFieldData().size() / 3;
404 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col,
false);
408 CHKERR iNtegrate(*(dataAtSpringPts->faceRowData), col_data);
411 CHKERR aSsemble(*(dataAtSpringPts->faceRowData), col_data);
423 const int *row_indices = &*row_data.
getIndices().data().begin();
425 const int *col_indices = &*col_data.
getIndices().data().begin();
427 auto &data = *dataAtSpringPts;
429 rowIndices.resize(row_nb_dofs,
false);
431 row_indices = &rowIndices[0];
433 VectorDofs::iterator dit = dofs.begin();
434 for (
int ii = 0; dit != dofs.end(); ++dit, ++ii) {
435 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
436 data.forcesOnlyOnEntitiesRow.end()) {
441 Mat
B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
442 : getFEMethod()->snes_B;
445 &*NN.data().begin(), ADD_VALUES);
460 auto get_tensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
462 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 1,
c + 0),
463 &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 2,
c + 0), &
m(r + 2,
c + 1),
470 const double normal_stiffness = dataAtSpringPts->springStiffnessNormal;
471 const double tangent_stiffness = dataAtSpringPts->springStiffnessTangent;
474 auto t_solution_at_gauss_point =
475 getFTensor1FromMat<3>(*dataAtSpringPts->xAtPts);
476 auto t_init_solution_at_gauss_point =
477 getFTensor1FromMat<3>(*dataAtSpringPts->xInitAtPts);
480 auto t_normal_1 = getFTensor1FromMat<3>(dataAtSpringPts->normalVector);
481 auto t_w = getFTensor0IntegrationWeight();
483 auto t_inv_H = getFTensor2FromMat<3, 3>(*dataAtSpringPts->invHMat);
484 auto t_F = getFTensor2FromMat<3, 3>(*dataAtSpringPts->FMat);
485 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
487 t_normal(
i) = t_normal_1(
i);
489 const double normal_length = std::sqrt(t_normal(
k) * t_normal(
k));
490 t_normal(
i) = t_normal(
i) / normal_length;
491 t_normal_projection(
i,
j) = t_normal(
i) * t_normal(
j);
492 t_tangent_projection(
i,
j) = t_normal_projection(
i,
j);
493 t_tangent_projection(0, 0) -= 1;
494 t_tangent_projection(1, 1) -= 1;
495 t_tangent_projection(2, 2) -= 1;
497 t_displacement_at_gauss_point(
i) =
498 t_solution_at_gauss_point(
i) - t_init_solution_at_gauss_point(
i);
500 double w = 0.5 * t_w * normal_length;
505 for (; bbc != nb_base_fun_col; ++bbc) {
510 for (; bbr != nb_base_fun_row; ++bbr) {
512 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
514 t_assemble(
i,
j) -= w * t_row_base * t_col_base * normal_stiffness *
515 t_F(
k,
i) * t_normal_projection(
k,
j);
516 t_assemble(
i,
j) += w * t_row_base * t_col_base * tangent_stiffness *
517 t_F(
k,
i) * t_tangent_projection(
k,
j);
519 t_assemble(
i,
j) -= w * normal_stiffness * t_row_base * t_inv_H(
k,
i) *
520 t_col_diff_base(
k) * t_normal_projection(
j,
l) *
521 t_displacement_at_gauss_point(
l);
522 t_assemble(
i,
j) += w * tangent_stiffness * t_row_base * t_inv_H(
k,
i) *
523 t_col_diff_base(
k) * t_tangent_projection(
j,
l) *
524 t_displacement_at_gauss_point(
l);
532 ++t_solution_at_gauss_point;
533 ++t_init_solution_at_gauss_point;
549 const int *row_indices = &*row_data.
getIndices().data().begin();
551 const int *col_indices = &*col_data.
getIndices().data().begin();
553 auto &data = *dataAtSpringPts;
555 rowIndices.resize(row_nb_dofs,
false);
557 row_indices = &rowIndices[0];
559 VectorDofs::iterator dit = dofs.begin();
560 for (
int ii = 0; dit != dofs.end(); ++dit, ++ii) {
561 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
562 data.forcesOnlyOnEntitiesRow.end()) {
567 Mat
B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
568 : getFEMethod()->snes_B;
571 &*NN.data().begin(), ADD_VALUES);
588 nb_gauss_pts = row_data.
getN().size1();
593 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col,
false);
596 diagonal_block = (row_type == col_type) && (row_side == col_side);
598 if (col_type == MBVERTEX) {
599 dataAtSpringPts->faceRowData = &row_data;
600 CHKERR loopSideVolumes(sideFeName, *sideFe);
604 CHKERR iNtegrate(row_data, col_data);
607 CHKERR aSsemble(row_data, col_data);
623 auto get_tensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
625 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 1,
c + 0),
626 &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 2,
c + 0), &
m(r + 2,
c + 1),
630 auto make_vec_der = [&](
auto t_N,
auto t_1,
auto t_2) {
638 auto make_vec_der_2 = [&](
auto t_N,
auto t_1,
auto t_2) {
641 const double normal_norm = std::sqrt(t_normal(
i) * t_normal(
i));
649 t_n(
i,
j) / normal_norm - t_normal(
i) * t_n(
k,
j) * t_normal(
k) /
650 (normal_norm * normal_norm * normal_norm);
654 dataAtSpringPts->faceRowData =
nullptr;
655 CHKERR loopSideVolumes(sideFeName, *sideFe);
657 auto t_F = getFTensor2FromMat<3, 3>(*dataAtSpringPts->FMat);
659 auto t_w = getFTensor0IntegrationWeight();
660 auto t_1 = getFTensor1FromMat<3>(dataAtSpringPts->tangent1);
661 auto t_2 = getFTensor1FromMat<3>(dataAtSpringPts->tangent2);
663 const double normal_stiffness = dataAtSpringPts->springStiffnessNormal;
664 const double tangent_stiffness = dataAtSpringPts->springStiffnessTangent;
673 auto t_solution_at_gauss_point =
674 getFTensor1FromMat<3>(*dataAtSpringPts->xAtPts);
675 auto t_init_solution_at_gauss_point =
676 getFTensor1FromMat<3>(*dataAtSpringPts->xInitAtPts);
679 auto t_normal_1 = getFTensor1FromMat<3>(dataAtSpringPts->normalVector);
682 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
683 t_normal(
i) = t_normal_1(
i);
685 const double normal_length = std::sqrt(t_normal(
k) * t_normal(
k));
686 t_normal(
i) = t_normal(
i) / normal_length;
687 t_normal_projection(
i,
j) = t_normal(
i) * t_normal(
j);
688 t_tangent_projection(
i,
j) = t_normal_projection(
i,
j);
689 t_tangent_projection(0, 0) -= 1;
690 t_tangent_projection(1, 1) -= 1;
691 t_tangent_projection(2, 2) -= 1;
693 t_displacement_at_gauss_point(
i) =
694 t_solution_at_gauss_point(
i) - t_init_solution_at_gauss_point(
i);
695 double val = 0.5 * t_w;
700 for (; bbc != nb_base_fun_col; ++bbc) {
705 for (; bbr != nb_base_fun_row; ++bbr) {
707 auto t_d_n = make_vec_der(t_N, t_1, t_2);
709 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
711 auto t_d_n_2 = make_vec_der_2(t_N, t_1, t_2);
713 t_assemble(
i,
j) += val * normal_length * t_col_base_func *
714 normal_stiffness * t_row_base_func * t_F(
k,
i) *
715 t_normal_projection(
k,
j);
717 t_assemble(
i,
j) -= val * normal_length * t_col_base_func *
718 tangent_stiffness * t_row_base_func * t_F(
k,
i) *
719 t_tangent_projection(
k,
j);
722 val * (normal_stiffness - tangent_stiffness) * t_row_base_func *
724 (t_d_n(
l,
j) * (t_normal(
k) * t_displacement_at_gauss_point(
k)) +
725 normal_length * t_normal(
l) *
726 (t_d_n_2(
k,
j) * t_displacement_at_gauss_point(
k)));
728 t_assemble(
i,
j) -= val * (t_d_n(
l,
j) * t_normal(
l)) *
729 tangent_stiffness * t_row_base_func * t_F(
k,
i) *
730 t_displacement_at_gauss_point(
k);
742 ++t_solution_at_gauss_point;
743 ++t_init_solution_at_gauss_point;
761 auto get_tensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
763 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 1,
c + 0),
764 &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 2,
c + 0), &
m(r + 2,
c + 1),
768 const double normal_stiffness = dataAtSpringPts->springStiffnessNormal;
769 const double tangent_stiffness = dataAtSpringPts->springStiffnessTangent;
778 auto t_solution_at_gauss_point =
779 getFTensor1FromMat<3>(*dataAtSpringPts->xAtPts);
780 auto t_init_solution_at_gauss_point =
781 getFTensor1FromMat<3>(*dataAtSpringPts->xInitAtPts);
784 auto t_normal_1 = getFTensor1FromMat<3>(dataAtSpringPts->normalVector);
786 auto t_w = getFTensor0IntegrationWeight();
788 auto t_h = getFTensor2FromMat<3, 3>(*dataAtSpringPts->hMat);
789 auto t_inv_H = getFTensor2FromMat<3, 3>(*dataAtSpringPts->invHMat);
793 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
794 t_normal(
i) = t_normal_1(
i);
796 const double normal_length = std::sqrt(t_normal(
k) * t_normal(
k));
797 t_normal(
i) = t_normal(
i) / normal_length;
798 t_normal_projection(
i,
j) = t_normal(
i) * t_normal(
j);
799 t_tangent_projection(
i,
j) = t_normal_projection(
i,
j);
800 t_tangent_projection(0, 0) -= 1;
801 t_tangent_projection(1, 1) -= 1;
802 t_tangent_projection(2, 2) -= 1;
804 t_displacement_at_gauss_point(
i) =
805 t_solution_at_gauss_point(
i) - t_init_solution_at_gauss_point(
i);
807 double val = 0.5 * t_w * normal_length;
812 for (; bbc != nb_base_fun_col; ++bbc) {
817 for (; bbr != nb_base_fun_row; ++bbr) {
819 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
822 val * t_row_base * normal_stiffness * t_inv_H(
l,
j) *
823 t_col_diff_base(
m) * t_inv_H(
m,
i) * t_h(
k,
l) *
824 (t_normal_projection(
k, q) * t_displacement_at_gauss_point(q));
827 val * t_row_base * tangent_stiffness * t_inv_H(
l,
j) *
828 t_col_diff_base(
m) * t_inv_H(
m,
i) * t_h(
k,
l) *
829 (t_tangent_projection(
k, q) * t_displacement_at_gauss_point(q));
838 ++t_solution_at_gauss_point;
839 ++t_init_solution_at_gauss_point;
856 const int nb_integration_pts = getGaussPts().size2();
858 auto t_h = getFTensor2FromMat<3, 3>(*commonDataPtr->hMat);
859 auto t_H = getFTensor2FromMat<3, 3>(*commonDataPtr->HMat);
861 commonDataPtr->detHVec->resize(nb_integration_pts,
false);
862 commonDataPtr->invHMat->resize(9, nb_integration_pts,
false);
863 commonDataPtr->FMat->resize(9, nb_integration_pts,
false);
865 commonDataPtr->detHVec->clear();
866 commonDataPtr->invHMat->clear();
867 commonDataPtr->FMat->clear();
870 auto t_invH = getFTensor2FromMat<3, 3>(*commonDataPtr->invHMat);
871 auto t_F = getFTensor2FromMat<3, 3>(*commonDataPtr->FMat);
873 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
876 t_F(
i,
j) = t_h(
i,
k) * t_invH(
k,
j);
893 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
904 nF.resize(nbRows,
false);
908 nbIntegrationPts = getGaussPts().size2();
911 CHKERR iNtegrate(row_data);
914 CHKERR aSsemble(row_data);
922 CHKERR loopSideVolumes(sideFeName, *sideFe);
929 CHKERR dataAtPts->getBlockData(dAta);
932 auto t_w = getFTensor0IntegrationWeight();
939 const double normal_stiffness = dataAtPts->springStiffnessNormal;
940 const double tangent_stiffness = dataAtPts->springStiffnessTangent;
949 auto t_solution_at_gauss_point = getFTensor1FromMat<3>(*dataAtPts->xAtPts);
950 auto t_init_solution_at_gauss_point =
951 getFTensor1FromMat<3>(*dataAtPts->xInitAtPts);
954 auto t_normal_1 = getFTensor1FromMat<3>(dataAtPts->normalVector);
956 auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
959 for (
int gg = 0; gg != nbIntegrationPts; ++gg) {
960 t_normal(
i) = t_normal_1(
i);
962 const double normal_length = std::sqrt(t_normal(
k) * t_normal(
k));
963 t_normal(
i) = t_normal(
i) / normal_length;
964 t_normal_projection(
i,
j) = t_normal(
i) * t_normal(
j);
965 t_tangent_projection(
i,
j) = t_normal_projection(
i,
j);
966 t_tangent_projection(0, 0) -= 1;
967 t_tangent_projection(1, 1) -= 1;
968 t_tangent_projection(2, 2) -= 1;
970 t_displacement_at_gauss_point(
i) =
971 t_solution_at_gauss_point(
i) - t_init_solution_at_gauss_point(
i);
973 double w = t_w * 0.5 * normal_length;
982 for (
int rr = 0; rr != nb_dofs / 3; ++rr) {
984 t_nf(
i) -= w * t_base_func * normal_stiffness * t_F(
k,
i) *
985 t_normal_projection(
k,
j) * t_displacement_at_gauss_point(
j);
986 t_nf(
i) += w * t_base_func * tangent_stiffness * t_F(
k,
i) *
987 t_tangent_projection(
k,
j) * t_displacement_at_gauss_point(
j);
997 ++t_solution_at_gauss_point;
998 ++t_init_solution_at_gauss_point;
1013 const int *row_indices = &*row_data.
getIndices().data().begin();
1015 auto &data = *dataAtPts;
1017 rowIndices.resize(nbRows,
false);
1019 row_indices = &rowIndices[0];
1021 VectorDofs::iterator dit = dofs.begin();
1022 for (
int ii = 0; dit != dofs.end(); ++dit, ++ii) {
1023 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
1024 data.forcesOnlyOnEntitiesRow.end()) {
1025 rowIndices[ii] = -1;
1029 CHKERR VecSetOption(getSNESf(), VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
1042 PetscFunctionReturn(0);
1044 ngp = data.
getN().size1();
1048 if (type == MBVERTEX) {
1049 dataAtIntegrationPts->tangent1.resize(3, ngp,
false);
1050 dataAtIntegrationPts->tangent1.clear();
1052 dataAtIntegrationPts->tangent2.resize(3, ngp,
false);
1053 dataAtIntegrationPts->tangent2.clear();
1056 auto t_1 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent1);
1057 auto t_2 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent2);
1059 for (
unsigned int gg = 0; gg != ngp; ++gg) {
1063 for (
unsigned int dd = 0; dd != nb_dofs; ++dd) {
1064 t_1(
i) += t_dof(
i) * t_N(0);
1065 t_2(
i) += t_dof(
i) * t_N(1);
1083 ngp = data.
getN().size1();
1089 if (type == MBVERTEX) {
1090 dataAtIntegrationPts->normalVector.resize(3, ngp,
false);
1091 dataAtIntegrationPts->normalVector.clear();
1094 auto normal_original_slave =
1095 getFTensor1FromMat<3>(dataAtIntegrationPts->normalVector);
1097 auto tangent_1 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent1);
1098 auto tangent_2 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent2);
1100 for (
unsigned int gg = 0; gg != ngp; ++gg) {
1101 normal_original_slave(
i) =
1103 ++normal_original_slave;
1117 if (col_type != MBVERTEX)
1120 dataAtSpringPts->faceRowData = &row_data;
1121 CHKERR loopSideVolumes(sideFeName, *sideFe);
1129 const std::string mesh_nodals_positions) {
1140 mesh_nodals_positions);
1145 if (
bit->getName().compare(0, 9,
"SPRING_BC") == 0) {
1156 const std::string mesh_nodals_positions,
Range &spring_triangles) {
1166 mesh_nodals_positions);
1168 mesh_nodals_positions);
1170 mesh_nodals_positions);
1180 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr,
1181 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr,
1182 const std::string
field_name,
const std::string mesh_nodals_positions,
double stiffness_scale) {
1187 boost::shared_ptr<MetaSpringBC::DataAtIntegrationPtsSprings> commonDataPtr =
1188 boost::make_shared<MetaSpringBC::DataAtIntegrationPtsSprings>(m_field, stiffness_scale);
1189 CHKERR commonDataPtr->getParameters();
1198 for (
auto &sitSpring : commonDataPtr->mapSpring) {
1200 fe_spring_lhs_ptr->getOpPtrVector().push_back(
1202 fe_spring_lhs_ptr->getOpPtrVector().push_back(
1204 fe_spring_lhs_ptr->getOpPtrVector().push_back(
1207 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1209 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1211 commonDataPtr->xInitAtPts));
1212 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1214 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1216 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1225 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr_dx,
1226 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr_dX,
1227 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr,
1228 boost::shared_ptr<MetaSpringBC::DataAtIntegrationPtsSprings>
1229 data_at_integration_pts,
1230 const std::string
field_name,
const std::string mesh_nodals_positions,
1231 std::string side_fe_name) {
1236 CHKERR data_at_integration_pts->getParameters();
1247 for (
auto &sitSpring : data_at_integration_pts->mapSpring) {
1249 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideRhs =
1250 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(m_field);
1252 feMatSideRhs->getOpPtrVector().push_back(
1254 mesh_nodals_positions, data_at_integration_pts->HMat));
1255 feMatSideRhs->getOpPtrVector().push_back(
1260 mesh_nodals_positions, data_at_integration_pts));
1262 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1264 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1267 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1269 data_at_integration_pts->xAtPts));
1270 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1272 mesh_nodals_positions, data_at_integration_pts->xInitAtPts));
1274 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1276 feMatSideRhs, side_fe_name, sitSpring.second));
1278 fe_spring_lhs_ptr_dx->getOpPtrVector().push_back(
1280 data_at_integration_pts->xAtPts));
1281 fe_spring_lhs_ptr_dx->getOpPtrVector().push_back(
1283 mesh_nodals_positions, data_at_integration_pts->xInitAtPts));
1284 fe_spring_lhs_ptr_dx->getOpPtrVector().push_back(
1286 fe_spring_lhs_ptr_dx->getOpPtrVector().push_back(
1289 fe_spring_lhs_ptr_dx->getOpPtrVector().push_back(
1291 mesh_nodals_positions));
1293 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideLhs_dx =
1294 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(m_field);
1296 feMatSideLhs_dx->getOpPtrVector().push_back(
1298 mesh_nodals_positions, data_at_integration_pts->HMat));
1299 feMatSideLhs_dx->getOpPtrVector().push_back(
1304 mesh_nodals_positions, data_at_integration_pts));
1306 feMatSideLhs_dx->getOpPtrVector().push_back(
1308 mesh_nodals_positions,
field_name, data_at_integration_pts));
1310 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1312 data_at_integration_pts->xAtPts));
1313 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1315 mesh_nodals_positions, data_at_integration_pts->xInitAtPts));
1316 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1318 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1321 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1323 data_at_integration_pts,
1324 feMatSideLhs_dx, side_fe_name));
1326 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideLhs_dX =
1327 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(m_field);
1329 feMatSideLhs_dX->getOpPtrVector().push_back(
1331 mesh_nodals_positions, data_at_integration_pts->HMat));
1332 feMatSideLhs_dX->getOpPtrVector().push_back(
1337 mesh_nodals_positions, data_at_integration_pts));
1339 feMatSideLhs_dX->getOpPtrVector().push_back(
1341 mesh_nodals_positions,
1342 data_at_integration_pts));
1344 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1346 mesh_nodals_positions, mesh_nodals_positions,
1347 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.
FTensor::Index< 'm', SPACE_DIM > m
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_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
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_row(const std::string &fe_name, const std::string name_row)=0
set field row 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
const double c
speed of light (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
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
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.
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
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.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
constexpr auto field_name
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.