12using namespace boost::numeric;
19 force = normal *
p / norm_2(normal);
28 boost::ptr_vector<MethodForForceScaling> &methods_op,
bool ho_geometry)
31 F(
_F), dAta(data),
methodsOp(methods_op), hoGeometry(ho_geometry) {}
40 EntityHandle ent = getNumeredEntFiniteElementPtr()->getEnt();
41 if (dAta.tRis.find(ent) == dAta.tRis.end())
45 int nb_row_dofs = data.
getIndices().size() / rank;
50 EntityType fe_type = getNumeredEntFiniteElementPtr()->getEntType();
52 for (
unsigned int gg = 0; gg != data.
getN().size1(); ++gg) {
55 double val = getGaussPts()(2, gg);
56 if (hoGeometry || fe_type == MBQUAD) {
57 val *= cblas_dnrm2(3, &getNormalsAtGaussPts()(gg, 0), 1);
64 for (
int rr = 0; rr != rank; ++rr) {
67 force = dAta.data.data.value3;
69 force = dAta.data.data.value4;
71 force = dAta.data.data.value5;
74 "data inconsistency");
76 force *= dAta.data.data.value1;
77 cblas_daxpy(nb_row_dofs, val * force, &data.
getN()(gg, 0), 1, &Nf[rr],
99 boost::ptr_vector<MethodForForceScaling> &methods_op,
100 boost::shared_ptr<MethodForAnalyticalForce> &analytical_force_op,
101 const bool ho_geometry)
113 const EntityHandle ent = getNumeredEntFiniteElementPtr()->getEnt();
114 if (tRis.find(ent) == tRis.end())
117 const int rank = data.
getFieldDofs()[0]->getNbOfCoeffs();
118 const int nb_row_dofs = data.
getIndices().size() / rank;
123 EntityType fe_type = getNumeredEntFiniteElementPtr()->getEntType();
129 for (
unsigned int gg = 0; gg != data.
getN().size1(); ++gg) {
132 double val = getGaussPts()(2, gg);
133 val *= cblas_dnrm2(3, &getNormalsAtGaussPts()(gg, 0), 1);
134 if (fe_type == MBTRI)
136 for (
int dd = 0; dd != 3; ++dd) {
137 coords[dd] = getCoordsAtGaussPts()(gg, dd);
138 normal[dd] = getNormalsAtGaussPts()(gg, dd);
143 for (
int rr = 0; rr != 3; ++rr) {
144 cblas_daxpy(nb_row_dofs, val * force[rr], &data.
getN()(gg, 0), 1,
169 boost::ptr_vector<MethodForForceScaling> &methods_op,
bool ho_geometry)
172 F(
_F), dAta(data),
methodsOp(methods_op), hoGeometry(ho_geometry) {}
182 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
183 EntityType fe_type = getNumeredEntFiniteElementPtr()->getEntType();
184 if (dAta.tRis.find(fe_ent) == dAta.tRis.end())
188 int nb_row_dofs = data.
getIndices().size() / rank;
193 for (
unsigned int gg = 0; gg != data.
getN().size1(); ++gg) {
195 double val = getGaussPts()(2, gg);
196 if (fe_type == MBTRI)
199 for (
int rr = 0; rr != rank; ++rr) {
202 if (hoGeometry || fe_type == MBQUAD)
203 force = dAta.data.data.value1 * getNormalsAtGaussPts()(gg, rr);
205 force = dAta.data.data.value1 * getNormal()[rr];
207 cblas_daxpy(nb_row_dofs, val * force, &data.
getN()(gg, 0), 1, &Nf[rr],
231 PetscFunctionReturn(0);
233 ngp = data.
getN().size1();
237 if (type == MBVERTEX) {
238 dataAtIntegrationPts->tangent1.resize(3, ngp,
false);
239 dataAtIntegrationPts->tangent1.clear();
241 dataAtIntegrationPts->tangent2.resize(3, ngp,
false);
242 dataAtIntegrationPts->tangent2.clear();
245 auto t_1 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent1);
246 auto t_2 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent2);
248 for (
unsigned int gg = 0; gg != ngp; ++gg) {
252 for (
unsigned int dd = 0; dd != nb_dofs; ++dd) {
253 t_1(
i) += t_dof(
i) * t_N(0);
254 t_2(
i) += t_dof(
i) * t_N(1);
271 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
276 const int row_nb_dofs = row_data.
getIndices().size();
279 const int col_nb_dofs = col_data.
getIndices().size();
282 const int nb_gauss_pts = row_data.
getN().size1();
284 int nb_base_fun_row = row_data.
getFieldData().size() / 3;
285 int nb_base_fun_col = col_data.
getFieldData().size() / 3;
287 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col,
false);
294 auto get_tensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
296 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 1,
c + 0),
297 &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 2,
c + 0), &
m(r + 2,
c + 1),
301 auto make_vec_der = [](
auto t_N,
auto t_1,
auto t_2) {
313 if (
auto arc_dof = dataAtIntegrationPts->arcLengthDof.lock()) {
314 lambda = arc_dof->getFieldData();
317 auto t_w = getFTensor0IntegrationWeight();
318 auto t_1 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent1);
319 auto t_2 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent2);
320 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
322 double val = 0.5 * t_w *
lambda;
327 for (; bbc != nb_base_fun_col; ++bbc) {
331 for (; bbr != nb_base_fun_row; ++bbr) {
333 auto t_d_n = make_vec_der(t_N, t_1, t_2);
335 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
338 t_assemble(
i,
k) += val * dAta.data.data.value1 * t_base * t_d_n(
i,
k);
350 const int *row_indices = &*row_data.
getIndices().data().begin();
352 const int *col_indices = &*col_data.
getIndices().data().begin();
354 Mat
B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
355 : getFEMethod()->snes_B;
358 &*NN.data().begin(), ADD_VALUES);
369 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
374 const int row_nb_dofs = row_data.
getIndices().size();
377 const int col_nb_dofs = col_data.
getIndices().size();
380 const int nb_gauss_pts = row_data.
getN().size1();
382 int nb_base_fun_row = row_data.
getFieldData().size() / 3;
383 int nb_base_fun_col = col_data.
getFieldData().size() / 3;
385 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col,
false);
392 auto get_tensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
394 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 1,
c + 0),
395 &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 2,
c + 0), &
m(r + 2,
c + 1),
399 auto make_vec_der = [](
auto t_N,
auto t_1,
auto t_2) {
411 if (
auto arc_dof = dataAtIntegrationPts->arcLengthDof.lock()) {
412 lambda = arc_dof->getFieldData();
415 auto t_w = getFTensor0IntegrationWeight();
416 auto t_1 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent1);
417 auto t_2 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent2);
418 auto t_normal = getFTensor1NormalsAtGaussPts();
419 FTensor1 t_force(dAta.data.data.value3, dAta.data.data.value4,
420 dAta.data.data.value5);
421 t_force(
i) *= dAta.
data.data.value1;
423 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
425 double nrm2 = std::sqrt(t_normal(
i) * t_normal(
i));
427 double val = 0.5 * t_w *
lambda / nrm2;
432 for (; bbc != nb_base_fun_col; ++bbc) {
436 for (; bbr != nb_base_fun_row; ++bbr) {
438 auto t_d_n = make_vec_der(t_N, t_1, t_2);
440 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
444 val * t_force(
i) * t_base * t_d_n(
j,
k) * t_normal(
j);
457 const int *row_indices = &*row_data.
getIndices().data().begin();
459 const int *col_indices = &*col_data.
getIndices().data().begin();
461 Mat
B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
462 : getFEMethod()->snes_B;
465 &*NN.data().begin(), ADD_VALUES);
479 const int nb_integration_pts = getGaussPts().size2();
480 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hMat);
481 auto t_H = getFTensor2FromMat<3, 3>(dataAtPts->HMat);
483 dataAtPts->detHVec.resize(nb_integration_pts,
false);
484 dataAtPts->invHMat.resize(9, nb_integration_pts,
false);
485 dataAtPts->FMat.resize(9, nb_integration_pts,
false);
488 auto t_invH = getFTensor2FromMat<3, 3>(dataAtPts->invHMat);
489 auto t_F = getFTensor2FromMat<3, 3>(dataAtPts->FMat);
491 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
494 t_F(
i,
j) = t_h(
i,
k) * t_invH(
k,
j);
511 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
522 nF.resize(nbRows,
false);
526 nbIntegrationPts = getGaussPts().size2();
529 CHKERR iNtegrate(row_data);
532 CHKERR aSsemble(row_data);
543 &
v(r + 0), &
v(r + 1), &
v(r + 2));
549 CHKERR loopSideVolumes(sideFeName, *sideFe);
552 auto t_w = getFTensor0IntegrationWeight();
554 auto t_normal = getFTensor1NormalsAtGaussPts();
556 auto t_F = getFTensor2FromMat<3, 3>(dataAtPts->FMat);
559 for (
int gg = 0; gg != nbIntegrationPts; ++gg) {
563 double a = -0.5 * t_w * dAta.data.data.value1;
564 auto t_nf = get_tensor1(nF, 0);
567 for (; rr != nbRows / 3; ++rr) {
569 t_nf(
i) +=
a * t_base * t_F(
j,
i) * t_normal(
j);
588 const int *row_indices = &*row_data.
getIndices().data().begin();
590 auto &data = *dataAtPts;
592 rowIndices.resize(nbRows,
false);
594 row_indices = &rowIndices[0];
596 VectorDofs::iterator dit = dofs.begin();
597 for (
int ii = 0; dit != dofs.end(); ++dit, ++ii) {
598 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
599 data.forcesOnlyOnEntitiesRow.end()) {
610 auto vec_assemble = [&](Vec my_f) {
612 CHKERR VecSetOption(my_f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
630 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
641 nF.resize(nbRows,
false);
645 nbIntegrationPts = getGaussPts().size2();
648 CHKERR iNtegrate(row_data);
651 CHKERR aSsemble(row_data);
663 &
v(r + 0), &
v(r + 1), &
v(r + 2));
669 CHKERR loopSideVolumes(sideFeName, *sideFe);
672 auto t_w = getFTensor0IntegrationWeight();
674 auto t_normal = getFTensor1NormalsAtGaussPts();
676 auto t_F = getFTensor2FromMat<3, 3>(dataAtPts->FMat);
678 FTensor1 t_force(dAta.data.data.value3, dAta.data.data.value4,
679 dAta.data.data.value5);
680 t_force(
i) *= dAta.
data.data.value1;
682 for (
int gg = 0; gg != nbIntegrationPts; ++gg) {
686 double nrm2 = std::sqrt(t_normal(
i) * t_normal(
i));
688 double a = -0.5 * t_w * nrm2;
689 auto t_nf = get_tensor1(nF, 0);
692 for (; rr != nbRows / 3; ++rr) {
694 t_nf(
i) +=
a * t_base * t_F(
j,
i) * t_force(
j);
714 const int *row_indices = &*row_data.
getIndices().data().begin();
716 auto &data = *dataAtPts;
718 rowIndices.resize(nbRows,
false);
720 row_indices = &rowIndices[0];
722 VectorDofs::iterator dit = dofs.begin();
723 for (
int ii = 0; dit != dofs.end(); ++dit, ++ii) {
724 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
725 data.forcesOnlyOnEntitiesRow.end()) {
736 auto vec_assemble = [&](Vec my_f) {
738 CHKERR VecSetOption(my_f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
757 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
768 nb_gauss_pts = row_data.
getN().size1();
773 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col,
false);
776 diagonal_block = (row_type == col_type) && (row_side == col_side);
778 if (col_type == MBVERTEX) {
779 dataAtPts->faceRowData = &row_data;
780 CHKERR loopSideVolumes(sideFeName, *sideFe);
784 CHKERR iNtegrate(row_data, col_data);
787 CHKERR aSsemble(row_data, col_data);
802 auto get_tensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
804 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 1,
c + 0),
805 &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 2,
c + 0), &
m(r + 2,
c + 1),
809 auto make_vec_der = [](
auto t_N,
auto t_1,
auto t_2) {
820 dataAtPts->faceRowData =
nullptr;
821 CHKERR loopSideVolumes(sideFeName, *sideFe);
823 auto t_F = getFTensor2FromMat<3, 3>(dataAtPts->FMat);
826 if (
auto arc_dof = dataAtPts->arcLengthDof.lock()) {
827 lambda = arc_dof->getFieldData();
830 auto t_w = getFTensor0IntegrationWeight();
831 auto t_1 = getFTensor1FromMat<3>(dataAtPts->tangent1);
832 auto t_2 = getFTensor1FromMat<3>(dataAtPts->tangent2);
834 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
836 double val = 0.5 * t_w *
lambda;
841 for (; bbc != nb_base_fun_col; ++bbc) {
846 for (; bbr != nb_base_fun_row; ++bbr) {
848 auto t_d_n = make_vec_der(t_N, t_1, t_2);
850 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
854 val * dAta.data.data.value1 * t_base * t_F(
j,
i) * t_d_n(
j,
k);
875 const int *row_indices = &*row_data.
getIndices().data().begin();
877 const int *col_indices = &*col_data.
getIndices().data().begin();
879 auto &data = *dataAtPts;
881 rowIndices.resize(row_nb_dofs,
false);
883 row_indices = &rowIndices[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()) {
893 Mat
B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
894 : getFEMethod()->snes_B;
897 &*NN.data().begin(), ADD_VALUES);
910 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
921 nb_gauss_pts = row_data.
getN().size1();
926 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col,
false);
929 diagonal_block = (row_type == col_type) && (row_side == col_side);
931 if (col_type == MBVERTEX) {
932 dataAtPts->faceRowData = &row_data;
933 CHKERR loopSideVolumes(sideFeName, *sideFe);
937 CHKERR iNtegrate(row_data, col_data);
940 CHKERR aSsemble(row_data, col_data);
956 auto get_tensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
958 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 1,
c + 0),
959 &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 2,
c + 0), &
m(r + 2,
c + 1),
963 auto make_vec_der = [](
auto t_N,
auto t_1,
auto t_2) {
974 dataAtPts->faceRowData =
nullptr;
975 CHKERR loopSideVolumes(sideFeName, *sideFe);
977 auto t_F = getFTensor2FromMat<3, 3>(dataAtPts->FMat);
980 if (
auto arc_dof = dataAtPts->arcLengthDof.lock()) {
981 lambda = arc_dof->getFieldData();
984 auto t_w = getFTensor0IntegrationWeight();
985 auto t_1 = getFTensor1FromMat<3>(dataAtPts->tangent1);
986 auto t_2 = getFTensor1FromMat<3>(dataAtPts->tangent2);
988 auto t_normal = getFTensor1NormalsAtGaussPts();
989 FTensor1 t_force(dAta.data.data.value3, dAta.data.data.value4,
990 dAta.data.data.value5);
991 t_force(
i) *= dAta.
data.data.value1;
992 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
994 double nrm2 = std::sqrt(t_normal(
i) * t_normal(
i));
996 double val = 0.5 * t_w *
lambda / nrm2;
1001 for (; bbc != nb_base_fun_col; ++bbc) {
1006 for (; bbr != nb_base_fun_row; ++bbr) {
1008 auto t_d_n = make_vec_der(t_N, t_1, t_2);
1010 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
1014 val * t_base * t_force(
j) * t_F(
j,
i) * (t_d_n(
l,
k) * t_normal(
l));
1036 const int *row_indices = &*row_data.
getIndices().data().begin();
1038 const int *col_indices = &*col_data.
getIndices().data().begin();
1040 auto &data = *dataAtPts;
1042 rowIndices.resize(row_nb_dofs,
false);
1044 row_indices = &rowIndices[0];
1046 VectorDofs::iterator dit = dofs.begin();
1047 for (
int ii = 0; dit != dofs.end(); ++dit, ++ii) {
1048 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
1049 data.forcesOnlyOnEntitiesRow.end()) {
1050 rowIndices[ii] = -1;
1054 Mat
B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
1055 : getFEMethod()->snes_B;
1058 &*NN.data().begin(), ADD_VALUES);
1070 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1075 if (col_type != MBVERTEX)
1078 dataAtPts->faceRowData = &row_data;
1079 CHKERR loopSideVolumes(sideFeName, *sideFe);
1092 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1097 if (col_type != MBVERTEX)
1100 dataAtPts->faceRowData = &row_data;
1101 CHKERR loopSideVolumes(sideFeName, *sideFe);
1114 if (dataAtPts->faceRowData ==
nullptr)
1117 if (row_type != MBVERTEX)
1120 row_nb_dofs = dataAtPts->faceRowData->getIndices().size();
1127 nb_gauss_pts = dataAtPts->faceRowData->getN().size1();
1129 nb_base_fun_row = dataAtPts->faceRowData->getFieldData().size() / 3;
1132 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col,
false);
1136 CHKERR iNtegrate(*(dataAtPts->faceRowData), col_data);
1139 CHKERR aSsemble(*(dataAtPts->faceRowData), col_data);
1152 if (dataAtPts->faceRowData ==
nullptr)
1155 if (row_type != MBVERTEX)
1158 row_nb_dofs = dataAtPts->faceRowData->getIndices().size();
1165 nb_gauss_pts = dataAtPts->faceRowData->getN().size1();
1167 nb_base_fun_row = dataAtPts->faceRowData->getFieldData().size() / 3;
1170 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col,
false);
1174 CHKERR iNtegrate(*(dataAtPts->faceRowData), col_data);
1177 CHKERR aSsemble(*(dataAtPts->faceRowData), col_data);
1192 auto get_tensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
1194 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 1,
c + 0),
1195 &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 2,
c + 0), &
m(r + 2,
c + 1),
1199 auto t_w = getFTensor0IntegrationWeight();
1201 auto t_normal = getFTensor1NormalsAtGaussPts();
1203 auto t_inv_H = getFTensor2FromMat<3, 3>(dataAtPts->invHMat);
1206 if (
auto arc_dof = dataAtPts->arcLengthDof.lock()) {
1207 lambda = arc_dof->getFieldData();
1210 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1212 double a = -0.5 * t_w * dAta.data.data.value1 *
lambda;
1217 for (; bbc != nb_base_fun_col; ++bbc) {
1222 for (; bbr != nb_base_fun_row; ++bbr) {
1224 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
1228 a * t_row_base * t_inv_H(
k,
i) * t_col_diff_base(
k) * t_normal(
j);
1252 auto get_tensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
1254 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 1,
c + 0),
1255 &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 2,
c + 0), &
m(r + 2,
c + 1),
1259 auto t_w = getFTensor0IntegrationWeight();
1261 auto t_normal = getFTensor1NormalsAtGaussPts();
1263 auto t_inv_H = getFTensor2FromMat<3, 3>(dataAtPts->invHMat);
1266 if (
auto arc_dof = dataAtPts->arcLengthDof.lock()) {
1267 lambda = arc_dof->getFieldData();
1270 FTensor1 t_force(dAta.data.data.value3, dAta.data.data.value4,
1271 dAta.data.data.value5);
1272 t_force(
i) *= dAta.
data.data.value1;
1273 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1275 double nrm2 = std::sqrt(t_normal(
i) * t_normal(
i));
1277 double a = -0.5 * t_w * nrm2 *
lambda;
1282 for (; bbc != nb_base_fun_col; ++bbc) {
1287 for (; bbr != nb_base_fun_row; ++bbr) {
1289 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
1293 a * t_row_base * t_inv_H(
k,
i) * t_col_diff_base(
k) * t_force(
j);
1314 const int *row_indices = &*row_data.
getIndices().data().begin();
1316 const int *col_indices = &*col_data.
getIndices().data().begin();
1318 auto &data = *dataAtPts;
1319 rowIndices.resize(row_nb_dofs,
false);
1321 row_indices = &rowIndices[0];
1323 VectorDofs::iterator dit = dofs.begin();
1324 for (
int ii = 0; dit != dofs.end(); ++dit, ++ii) {
1325 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
1326 data.forcesOnlyOnEntitiesRow.end()) {
1327 rowIndices[ii] = -1;
1331 Mat
B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
1332 : getFEMethod()->snes_B;
1335 &*NN.data().begin(), ADD_VALUES);
1347 const int *row_indices = &*row_data.
getIndices().data().begin();
1349 const int *col_indices = &*col_data.
getIndices().data().begin();
1351 auto &data = *dataAtPts;
1353 rowIndices.resize(row_nb_dofs,
false);
1355 row_indices = &rowIndices[0];
1357 VectorDofs::iterator dit = dofs.begin();
1358 for (
int ii = 0; dit != dofs.end(); ++dit, ++ii) {
1359 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
1360 data.forcesOnlyOnEntitiesRow.end()) {
1361 rowIndices[ii] = -1;
1365 Mat
B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
1366 : getFEMethod()->snes_B;
1370 &*NN.data().begin(), ADD_VALUES);
1387 auto get_tensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
1389 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 1,
c + 0),
1390 &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 2,
c + 0), &
m(r + 2,
c + 1),
1394 auto t_w = getFTensor0IntegrationWeight();
1396 auto t_normal = getFTensor1NormalsAtGaussPts();
1398 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hMat);
1399 auto t_inv_H = getFTensor2FromMat<3, 3>(dataAtPts->invHMat);
1404 if (
auto arc_dof = dataAtPts->arcLengthDof.lock()) {
1405 lambda = arc_dof->getFieldData();
1408 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1410 double a = -0.5 * t_w * dAta.data.data.value1 *
lambda;
1415 for (; bbc != nb_base_fun_col; ++bbc) {
1420 for (; bbr != nb_base_fun_row; ++bbr) {
1422 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
1426 t_assemble(
i,
j) -=
a * t_row_base * t_inv_H(
l,
j) *
1427 t_col_diff_base(
m) * t_inv_H(
m,
i) * t_h(
k,
l) *
1455 auto get_tensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
1457 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 1,
c + 0),
1458 &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 2,
c + 0), &
m(r + 2,
c + 1),
1462 auto t_w = getFTensor0IntegrationWeight();
1464 auto t_normal = getFTensor1NormalsAtGaussPts();
1466 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hMat);
1467 auto t_inv_H = getFTensor2FromMat<3, 3>(dataAtPts->invHMat);
1472 if (
auto arc_dof = dataAtPts->arcLengthDof.lock()) {
1473 lambda = arc_dof->getFieldData();
1475 FTensor1 t_force(dAta.data.data.value3, dAta.data.data.value4,
1476 dAta.data.data.value5);
1477 t_force(
i) *= dAta.
data.data.value1;
1478 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1480 double nrm2 = std::sqrt(t_normal(
i) * t_normal(
i));
1481 double a = -0.5 * t_w * nrm2 *
lambda;
1486 for (; bbc != nb_base_fun_col; ++bbc) {
1491 for (; bbr != nb_base_fun_row; ++bbr) {
1493 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
1497 t_assemble(
i,
j) -=
a * t_row_base * t_inv_H(
l,
j) *
1498 t_col_diff_base(
m) * t_inv_H(
m,
i) * t_h(
k,
l) *
1516 boost::ptr_vector<MethodForForceScaling> &methods_op,
bool ho_geometry)
1519 F(
_F), dAta(data),
methodsOp(methods_op), hoGeometry(ho_geometry) {}
1528 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1534 int nb_row_dofs = data.
getIndices().size() / rank;
1539 EntityType fe_type = getNumeredEntFiniteElementPtr()->getEntType();
1541 for (
unsigned int gg = 0; gg != data.
getN().size1(); ++gg) {
1543 double val = getGaussPts()(2, gg);
1545 if (hoGeometry || fe_type == MBQUAD) {
1546 double area = cblas_dnrm2(3, &getNormalsAtGaussPts()(gg, 0), 1);
1547 if (fe_type == MBTRI)
1549 flux = dAta.data.data.value1 * area;
1551 flux = dAta.data.data.value1 * getArea();
1554 cblas_daxpy(nb_row_dofs, val * flux, &data.
getN()(gg, 0), 1,
1555 &*Nf.data().begin(), 1);
1560 auto get_f = [&]() {
1561 if (
F == PETSC_NULL)
1582 &cubit_meshset_ptr);
1583 std::vector<double> mydata;
1586 for (
unsigned int ii = 0; ii != mydata.size(); ++ii) {
1587 force[ii] = mydata[ii];
1590 if (force.empty()) {
1595 const string name =
"Force";
1596 strncpy(
mapForce[ms_id].data.data.name, name.c_str(),
1597 name.size() > 5 ? 5 : name.size());
1598 double magnitude = std::sqrt(force[0] * force[0] + force[1] * force[1] +
1599 force[2] * force[2]);
1600 mapForce[ms_id].data.data.value1 = -magnitude;
1601 mapForce[ms_id].data.data.value2 = 0;
1603 force[0] / magnitude;
1605 force[1] / magnitude;
1607 force[2] / magnitude;
1608 mapForce[ms_id].data.data.value6 = 0;
1609 mapForce[ms_id].data.data.value7 = 0;
1610 mapForce[ms_id].data.data.value8 = 0;
1611 mapForce[ms_id].data.data.zero[0] = 0;
1612 mapForce[ms_id].data.data.zero[1] = 0;
1613 mapForce[ms_id].data.data.zero[2] = 0;
1614 mapForce[ms_id].data.data.zero2 = 0;
1634 const std::string x_field,
const std::string X_field,
1635 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1636 std::string side_fe_name, Vec
F, Mat aij,
int ms_id,
bool ho_geometry,
1637 bool block_set,
bool ignore_material_force) {
1645 &cubit_meshset_ptr);
1646 std::vector<double> mydata;
1649 for (
unsigned int ii = 0; ii != mydata.size(); ++ii) {
1650 force[ii] = mydata[ii];
1653 if (force.empty()) {
1658 const string name =
"Force";
1659 strncpy(
mapForce[ms_id].data.data.name, name.c_str(),
1660 name.size() > 5 ? 5 : name.size());
1661 double magnitude = std::sqrt(force[0] * force[0] + force[1] * force[1] +
1662 force[2] * force[2]);
1663 mapForce[ms_id].data.data.value1 = -magnitude;
1664 mapForce[ms_id].data.data.value2 = 0;
1666 force[0] / magnitude;
1668 force[1] / magnitude;
1670 force[2] / magnitude;
1671 mapForce[ms_id].data.data.value6 = 0;
1672 mapForce[ms_id].data.data.value7 = 0;
1673 mapForce[ms_id].data.data.value8 = 0;
1674 mapForce[ms_id].data.data.zero[0] = 0;
1675 mapForce[ms_id].data.data.zero[1] = 0;
1676 mapForce[ms_id].data.data.zero[2] = 0;
1677 mapForce[ms_id].data.data.zero2 = 0;
1693 x_field, X_field, data_at_pts, aij,
mapForce[ms_id], ho_geometry));
1696 if (!ignore_material_force) {
1698 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideRhs =
1699 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(
mField);
1701 feMatSideRhs->getOpPtrVector().push_back(
1703 data_at_pts->getHMatPtr()));
1704 feMatSideRhs->getOpPtrVector().push_back(
1706 x_field, data_at_pts->getSmallhMatPtr()));
1708 feMatSideRhs->getOpPtrVector().push_back(
1712 X_field, data_at_pts, feMatSideRhs, side_fe_name,
F,
mapForce[ms_id],
1718 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideLhs_dx =
1719 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(
mField);
1721 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideLhs_dX =
1722 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(
mField);
1724 feMatSideLhs_dx->getOpPtrVector().push_back(
1726 data_at_pts->getHMatPtr()));
1727 feMatSideLhs_dx->getOpPtrVector().push_back(
1729 x_field, data_at_pts->getSmallhMatPtr()));
1731 feMatSideLhs_dx->getOpPtrVector().push_back(
1733 feMatSideLhs_dx->getOpPtrVector().push_back(
1735 X_field, x_field, data_at_pts, aij,
mapForce[ms_id], ho_geometry));
1739 X_field, x_field, data_at_pts, feMatSideLhs_dx, side_fe_name, aij,
1742 feMatSideLhs_dX->getOpPtrVector().push_back(
1744 data_at_pts->getHMatPtr()));
1745 feMatSideLhs_dX->getOpPtrVector().push_back(
1747 x_field, data_at_pts->getSmallhMatPtr()));
1748 feMatSideLhs_dX->getOpPtrVector().push_back(
1750 feMatSideLhs_dX->getOpPtrVector().push_back(
1752 X_field, X_field, data_at_pts, aij,
mapForce[ms_id], ho_geometry));
1757 X_field, X_field, data_at_pts, feMatSideLhs_dX, side_fe_name, aij,
1775 &cubit_meshset_ptr);
1776 std::vector<double> mydata;
1779 for (
unsigned int ii = 0; ii != mydata.size(); ++ii) {
1780 pressure[ii] = mydata[ii];
1782 if (pressure.empty()) {
1785 const string name =
"Pressure";
1786 strncpy(
mapPressure[ms_id].data.data.name, name.c_str(),
1787 name.size() > 8 ? 8 : name.size());
1790 mapPressure[ms_id].data.data.value1 = pressure[0];
1808 const std::string x_field,
const std::string X_field,
1809 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1810 std::string side_fe_name, Vec
F, Mat aij,
int ms_id,
bool ho_geometry,
1819 &cubit_meshset_ptr);
1820 std::vector<double> mydata;
1823 for (
unsigned int ii = 0; ii != mydata.size(); ++ii) {
1824 pressure[ii] = mydata[ii];
1826 if (pressure.empty()) {
1829 const string name =
"Pressure";
1830 strncpy(
mapPressure[ms_id].data.data.name, name.c_str(),
1831 name.size() > 8 ? 8 : name.size());
1834 mapPressure[ms_id].data.data.value1 = pressure[0];
1851 x_field, X_field, data_at_pts, aij,
mapPressure[ms_id], ho_geometry));
1856 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideRhs =
1857 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(
mField);
1859 feMatSideRhs->getOpPtrVector().push_back(
1861 data_at_pts->getHMatPtr()));
1862 feMatSideRhs->getOpPtrVector().push_back(
1864 data_at_pts->getSmallhMatPtr()));
1866 feMatSideRhs->getOpPtrVector().push_back(
1870 X_field, data_at_pts, feMatSideRhs, side_fe_name,
F,
mapPressure[ms_id],
1876 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideLhs_dx =
1877 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(
mField);
1879 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideLhs_dX =
1880 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(
mField);
1882 feMatSideLhs_dx->getOpPtrVector().push_back(
1884 data_at_pts->getHMatPtr()));
1885 feMatSideLhs_dx->getOpPtrVector().push_back(
1887 data_at_pts->getSmallhMatPtr()));
1889 feMatSideLhs_dx->getOpPtrVector().push_back(
1891 feMatSideLhs_dx->getOpPtrVector().push_back(
1893 X_field, x_field, data_at_pts, aij,
mapPressure[ms_id], ho_geometry));
1896 X_field, x_field, data_at_pts, feMatSideLhs_dx, side_fe_name, aij,
1899 feMatSideLhs_dX->getOpPtrVector().push_back(
1901 data_at_pts->getHMatPtr()));
1902 feMatSideLhs_dX->getOpPtrVector().push_back(
1904 data_at_pts->getSmallhMatPtr()));
1905 feMatSideLhs_dX->getOpPtrVector().push_back(
1907 feMatSideLhs_dX->getOpPtrVector().push_back(
1909 X_field, X_field, data_at_pts, aij,
mapPressure[ms_id], ho_geometry));
1913 X_field, X_field, data_at_pts, feMatSideLhs_dX, side_fe_name, aij,
1921 int ms_id,
bool ho_geometry) {
1928 std::vector<double> mydata;
1930 if (mydata.size() < 4)
1932 "Should be four block attributes but is %d", mydata.size());
1934 for (
unsigned int ii = 0; ii != 3; ++ii) {
1935 pressure_coeffs[ii] = mydata[ii];
1937 const double pressure_shift = mydata[3];
1942 boost::shared_ptr<MethodForAnalyticalForce> analytical_force_op(
1976 const std::string mesh_nodals_positions,
Range *intersect_ptr) {
1987 mesh_nodals_positions);
1994 CHKERR m_field.
get_moab().get_entities_by_dimension(it->meshset, 2, tris,
1997 tris = intersect(tris, *intersect_ptr);
2008 mesh_nodals_positions);
2014 CHKERR m_field.
get_moab().get_entities_by_dimension(it->meshset, 2, tris,
2017 tris = intersect(tris, *intersect_ptr);
2023 const string block_set_force_name(
"FORCE");
2026 if (it->getName().compare(0, block_set_force_name.length(),
2027 block_set_force_name) == 0) {
2029 CHKERR m_field.
get_moab().get_entities_by_dimension(it->meshset, 2, tris,
2032 tris = intersect(tris, *intersect_ptr);
2038 const string block_set_pressure_name(
"PRESSURE");
2040 if (it->getName().compare(0, block_set_pressure_name.length(),
2041 block_set_pressure_name) == 0) {
2043 CHKERR m_field.
get_moab().get_entities_by_dimension(it->meshset, 2, tris,
2046 tris = intersect(tris, *intersect_ptr);
2053 const string block_set_linear_pressure_name(
"LINEAR_PRESSURE");
2055 if (it->getName().compare(0, block_set_linear_pressure_name.length(),
2056 block_set_linear_pressure_name) == 0) {
2058 CHKERR m_field.
get_moab().get_entities_by_dimension(it->meshset, 2, tris,
2061 tris = intersect(tris, *intersect_ptr);
2071 boost::ptr_map<std::string, NeumannForcesSurface> &neumann_forces, Vec
F,
2072 const std::string
field_name,
const std::string mesh_nodals_positions) {
2074 bool ho_geometry = m_field.
check_field(mesh_nodals_positions);
2077 std::string fe_name =
"FORCE_FE";
2079 auto &nf = neumann_forces.at(fe_name);
2080 auto &
fe = nf.getLoopFe();
2090 const string block_set_force_name(
"FORCE");
2092 if (it->getName().compare(0, block_set_force_name.length(),
2093 block_set_force_name) == 0) {
2102 std::string fe_name =
"PRESSURE_FE";
2105 auto &nf = neumann_forces.at(fe_name);
2106 auto &
fe = nf.getLoopFe();
2118 const string block_set_pressure_name(
"PRESSURE");
2120 if (it->getName().compare(0, block_set_pressure_name.length(),
2121 block_set_pressure_name) == 0) {
2128 const string block_set_linear_pressure_name(
"LINEAR_PRESSURE");
2130 if (it->getName().compare(0, block_set_linear_pressure_name.length(),
2131 block_set_linear_pressure_name) == 0) {
2143 const std::string mesh_nodals_positions) {
2152 mesh_nodals_positions);
2158 CHKERR m_field.
get_moab().get_entities_by_dimension(it->meshset, 2, tris,
2168 boost::ptr_map<std::string, NeumannForcesSurface> &neumann_forces, Vec
F,
2169 const std::string
field_name,
const std::string mesh_nodals_positions) {
2171 bool ho_geometry = m_field.
check_field(mesh_nodals_positions);
2174 fe_name =
"FLUX_FE";
2177 auto &nf = neumann_forces.at(fe_name);
2178 auto &
fe = nf.getLoopFe();
#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 ...
@ MOFEM_DATA_INCONSISTENCY
#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_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
#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)
const double v
phase velocity of light in medium (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.
VectorBoundedArray< double, 3 > VectorDouble3
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
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
virtual moab::Interface & get_moab()=0
this struct keeps basic methods for moab meshset about material and boundary conditions
MoFEMErrorCode getAttributes(std::vector< double > &attributes) const
get Cubit block attributes
MoFEMErrorCode getBcDataStructure(CUBIT_BC_DATA_TYPE &data) const
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.
default operator for TRI element
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Interface for managing meshsets containing materials and boundary conditions.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MoFEMErrorCode getForce(const EntityHandle ent, const VectorDouble3 &coords, const VectorDouble3 &normal, VectorDouble3 &force)
const double pressureShift
const VectorDouble3 linearConstants
MyTriangleFE(MoFEM::Interface &m_field)
Operator for computing tangent vectors.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for flux element.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpNeumannFlux(const std::string field_name, Vec _F, bCPressure &data, boost::ptr_vector< MethodForForceScaling > &methods_op, bool ho_geometry)
Operator for force element.
OpNeumannForceAnalytical(const std::string field_name, Vec f, const Range tris, boost::ptr_vector< MethodForForceScaling > &methods_op, boost::shared_ptr< MethodForAnalyticalForce > &analytical_force_op, const bool ho_geometry=false)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for force element.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Integrate surface force (traction)
OpNeumannForce(const std::string field_name, Vec _F, bCForce &data, boost::ptr_vector< MethodForForceScaling > &methods_op, bool ho_geometry=false)
RHS-operator for pressure element (spatial configuration)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Integrate pressure.
OpNeumannPressure(const std::string field_name, Vec _F, bCPressure &data, boost::ptr_vector< MethodForForceScaling > &methods_op, bool ho_geometry=false)
LHS-operator for pressure element (spatial configuration)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Compute left-hand side.
LHS-operator for the pressure element (material configuration)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Compute part of the left-hand side.
LHS-operator for the pressure element (material configuration)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
RHS-operator for the pressure element (material configuration)
MoFEMErrorCode iNtegrate(EntData &row_data)
MoFEMErrorCode aSsemble(EntData &row_data)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &row_data)
Integrate pressure in the material configuration.
LHS-operator (material configuration) on the side volume.
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrates over a face contribution from a side volume.
LHS-operator (material configuration) on the side volume.
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrates over a face contribution from a side volume.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
LHS-operator for surface force element (spatial configuration)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Compute left-hand side.
LHS-operator for the surface force element (material configuration)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Compute part of the left-hand side.
LHS-operator for the surface force element (material configuration)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
RHS-operator for the surface force element (material configuration)
MoFEMErrorCode aSsemble(EntData &row_data)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &row_data)
Integrate surface force in the material configuration.
MoFEMErrorCode iNtegrate(EntData &row_data)
LHS-operator (material configuration) on the side volume.
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrates over a face contribution from a side volume.
LHS-operator (material configuration) on the side volume.
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrates over a face contribution from a side volume.
MoFEMErrorCode aSsemble(EntData &row_data, 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)
Finite element and operators to apply force/pressures applied to surfaces.
std::map< int, bCPressure > mapPressure
MoFEMErrorCode addPressureAle(const std::string field_name_1, const std::string field_name_2, boost::shared_ptr< DataAtIntegrationPts > data_at_pts, std::string side_fe_name, Vec F, Mat aij, int ms_id, bool ho_geometry=false, bool block_set=false)
Add operator to calculate pressure on element (in ALE)
MoFEMErrorCode addForceAle(const std::string field_name_1, const std::string field_name_2, boost::shared_ptr< DataAtIntegrationPts > data_at_pts, std::string side_fe_name, Vec F, Mat aij, int ms_id, bool ho_geometry=false, bool block_set=false, bool ignore_material_force=false)
Add operator to calculate forces on element (in ALE)
std::map< int, bCForce > mapForce
MoFEMErrorCode addFlux(const std::string field_name, Vec F, int ms_id, bool ho_geometry=false)
Add flux element operator (integration on face)
boost::ptr_vector< MethodForForceScaling > methodsOp
DEPRECATED MoFEMErrorCode addPreassure(const std::string field_name, Vec F, int ms_id, bool ho_geometry=false, bool block_set=false)
MoFEMErrorCode addForce(const std::string field_name, Vec F, int ms_id, bool ho_geometry=false, bool block_set=false)
Add operator to calculate forces on element.
MoFEMErrorCode addLinearPressure(const std::string field_name, Vec F, int ms_id, bool ho_geometry=false)
Add operator to calculate pressure on element.
MoFEMErrorCode addPressure(const std::string field_name, Vec F, int ms_id, bool ho_geometry=false, bool block_set=false)
Add operator to calculate pressure on element.
boost::ptr_vector< MethodForAnalyticalForce > analyticalForceOp
MoFEM::Interface & mField