24using namespace boost::numeric;
40 boost::ptr_vector<MethodForForceScaling> &methods_op,
bool ho_geometry)
43 F(
_F), dAta(data),
methodsOp(methods_op), hoGeometry(ho_geometry) {}
51 EntityHandle ent = getNumeredEntFiniteElementPtr()->getEnt();
52 if (dAta.tRis.find(ent) == dAta.tRis.end())
56 int nb_row_dofs = data.
getIndices().size() / rank;
61 EntityType fe_type = getNumeredEntFiniteElementPtr()->getEntType();
63 for (
unsigned int gg = 0; gg != data.
getN().size1(); ++gg) {
66 double val = getGaussPts()(2, gg);
67 if (hoGeometry || fe_type == MBQUAD) {
68 val *= cblas_dnrm2(3, &getNormalsAtGaussPts()(gg, 0), 1);
75 for (
int rr = 0; rr != rank; ++rr) {
78 force = dAta.data.data.value3;
80 force = dAta.data.data.value4;
82 force = dAta.data.data.value5;
85 "data inconsistency");
87 force *= dAta.data.data.value1;
88 cblas_daxpy(nb_row_dofs, val * force, &data.
getN()(gg, 0), 1, &Nf[rr],
110 boost::ptr_vector<MethodForForceScaling> &methods_op,
111 boost::shared_ptr<MethodForAnalyticalForce> &analytical_force_op,
112 const bool ho_geometry)
124 const EntityHandle ent = getNumeredEntFiniteElementPtr()->getEnt();
125 if (tRis.find(ent) == tRis.end())
128 const int rank = data.
getFieldDofs()[0]->getNbOfCoeffs();
129 const int nb_row_dofs = data.
getIndices().size() / rank;
134 EntityType fe_type = getNumeredEntFiniteElementPtr()->getEntType();
140 for (
unsigned int gg = 0; gg != data.
getN().size1(); ++gg) {
143 double val = getGaussPts()(2, gg);
144 val *= cblas_dnrm2(3, &getNormalsAtGaussPts()(gg, 0), 1);
145 if (fe_type == MBTRI)
147 for (
int dd = 0;
dd != 3; ++
dd) {
148 coords[
dd] = getCoordsAtGaussPts()(gg,
dd);
154 for (
int rr = 0; rr != 3; ++rr) {
155 cblas_daxpy(nb_row_dofs, val * force[rr], &data.
getN()(gg, 0), 1,
180 boost::ptr_vector<MethodForForceScaling> &methods_op,
bool ho_geometry)
183 F(
_F), dAta(data),
methodsOp(methods_op), hoGeometry(ho_geometry) {}
193 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
194 EntityType fe_type = getNumeredEntFiniteElementPtr()->getEntType();
195 if (dAta.tRis.find(fe_ent) == dAta.tRis.end())
199 int nb_row_dofs = data.
getIndices().size() / rank;
204 for (
unsigned int gg = 0; gg != data.
getN().size1(); ++gg) {
206 double val = getGaussPts()(2, gg);
207 if (fe_type == MBTRI)
210 for (
int rr = 0; rr != rank; ++rr) {
213 if (hoGeometry || fe_type == MBQUAD)
214 force = dAta.data.data.value1 * getNormalsAtGaussPts()(gg, rr);
216 force = dAta.data.data.value1 * getNormal()[rr];
218 cblas_daxpy(nb_row_dofs, val * force, &data.
getN()(gg, 0), 1, &Nf[rr],
241 PetscFunctionReturn(0);
243 ngp = data.
getN().size1();
247 if (
type == MBVERTEX) {
248 dataAtIntegrationPts->tangent1.resize(3, ngp,
false);
249 dataAtIntegrationPts->tangent1.clear();
251 dataAtIntegrationPts->tangent2.resize(3, ngp,
false);
252 dataAtIntegrationPts->tangent2.clear();
255 auto t_1 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent1);
256 auto t_2 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent2);
258 for (
unsigned int gg = 0; gg != ngp; ++gg) {
262 for (
unsigned int dd = 0;
dd != nb_dofs; ++
dd) {
263 t_1(
i) += t_dof(
i) * t_N(0);
264 t_2(
i) += t_dof(
i) * t_N(1);
281 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
286 const int row_nb_dofs = row_data.
getIndices().size();
289 const int col_nb_dofs = col_data.
getIndices().size();
292 const int nb_gauss_pts = row_data.
getN().size1();
294 int nb_base_fun_row = row_data.
getFieldData().size() / 3;
295 int nb_base_fun_col = col_data.
getFieldData().size() / 3;
297 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col,
false);
306 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2), &
m(
r + 1,
c + 0),
307 &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2), &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1),
311 auto make_vec_der = [](
auto t_N,
auto t_1,
auto t_2) {
323 if (
auto arc_dof = dataAtIntegrationPts->arcLengthDof.lock()) {
324 lambda = arc_dof->getFieldData();
327 auto t_w = getFTensor0IntegrationWeight();
328 auto t_1 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent1);
329 auto t_2 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent2);
330 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
332 double val = 0.5 * t_w *
lambda;
337 for (; bbc != nb_base_fun_col; ++bbc) {
341 for (; bbr != nb_base_fun_row; ++bbr) {
343 auto t_d_n = make_vec_der(t_N, t_1, t_2);
345 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
348 t_assemble(
i,
k) += val * dAta.data.data.value1 * t_base * t_d_n(
i,
k);
360 const int *row_indices = &*row_data.
getIndices().data().begin();
362 const int *col_indices = &*col_data.
getIndices().data().begin();
364 Mat
B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
365 : getFEMethod()->snes_B;
368 &*NN.data().begin(), ADD_VALUES);
382 const int nb_integration_pts = getGaussPts().size2();
383 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hMat);
384 auto t_H = getFTensor2FromMat<3, 3>(dataAtPts->HMat);
386 dataAtPts->detHVec.resize(nb_integration_pts,
false);
387 dataAtPts->invHMat.resize(9, nb_integration_pts,
false);
388 dataAtPts->FMat.resize(9, nb_integration_pts,
false);
391 auto t_invH = getFTensor2FromMat<3, 3>(dataAtPts->invHMat);
392 auto t_F = getFTensor2FromMat<3, 3>(dataAtPts->FMat);
394 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
397 t_F(
i,
j) = t_h(
i,
k) * t_invH(
k,
j);
414 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
425 nF.resize(nbRows,
false);
429 nbIntegrationPts = getGaussPts().size2();
432 CHKERR iNtegrate(row_data);
435 CHKERR aSsemble(row_data);
446 &
v(
r + 0), &
v(
r + 1), &
v(
r + 2));
452 CHKERR loopSideVolumes(sideFeName, *sideFe);
455 auto t_w = getFTensor0IntegrationWeight();
457 auto t_normal = getFTensor1NormalsAtGaussPts();
459 auto t_F = getFTensor2FromMat<3, 3>(dataAtPts->FMat);
462 for (
int gg = 0; gg != nbIntegrationPts; ++gg) {
466 double a = -0.5 * t_w * dAta.data.data.value1;
467 auto t_nf = get_tensor1(nF, 0);
470 for (; rr != nbRows / 3; ++rr) {
472 t_nf(
i) +=
a * t_base * t_F(
j,
i) * t_normal(
j);
491 const int *row_indices = &*row_data.
getIndices().data().begin();
493 auto &data = *dataAtPts;
494 if (!data.forcesOnlyOnEntitiesRow.empty()) {
495 rowIndices.resize(nbRows,
false);
497 row_indices = &rowIndices[0];
499 VectorDofs::iterator dit = dofs.begin();
500 for (
int ii = 0; dit != dofs.end(); ++dit, ++ii) {
501 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
502 data.forcesOnlyOnEntitiesRow.end()) {
514 auto vec_assemble = [&](
Vec my_f) {
516 CHKERR VecSetOption(my_f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
535 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
546 nb_gauss_pts = row_data.
getN().size1();
551 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col,
false);
554 diagonal_block = (row_type == col_type) && (row_side == col_side);
556 if (col_type == MBVERTEX) {
557 dataAtPts->faceRowData = &row_data;
558 CHKERR loopSideVolumes(sideFeName, *sideFe);
562 CHKERR iNtegrate(row_data, col_data);
565 CHKERR aSsemble(row_data, col_data);
582 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2), &
m(
r + 1,
c + 0),
583 &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2), &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1),
587 auto make_vec_der = [](
auto t_N,
auto t_1,
auto t_2) {
598 dataAtPts->faceRowData =
nullptr;
599 CHKERR loopSideVolumes(sideFeName, *sideFe);
601 auto t_F = getFTensor2FromMat<3, 3>(dataAtPts->FMat);
604 if (
auto arc_dof = dataAtPts->arcLengthDof.lock()) {
605 lambda = arc_dof->getFieldData();
608 auto t_w = getFTensor0IntegrationWeight();
609 auto t_1 = getFTensor1FromMat<3>(dataAtPts->tangent1);
610 auto t_2 = getFTensor1FromMat<3>(dataAtPts->tangent2);
612 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
614 double val = 0.5 * t_w *
lambda;
619 for (; bbc != nb_base_fun_col; ++bbc) {
624 for (; bbr != nb_base_fun_row; ++bbr) {
626 auto t_d_n = make_vec_der(t_N, t_1, t_2);
628 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
632 val * dAta.data.data.value1 * t_base * t_F(
j,
i) * t_d_n(
j,
k);
653 const int *row_indices = &*row_data.
getIndices().data().begin();
655 const int *col_indices = &*col_data.
getIndices().data().begin();
657 auto &data = *dataAtPts;
658 if (!data.forcesOnlyOnEntitiesRow.empty()) {
659 rowIndices.resize(row_nb_dofs,
false);
661 row_indices = &rowIndices[0];
663 VectorDofs::iterator dit = dofs.begin();
664 for (
int ii = 0; dit != dofs.end(); ++dit, ++ii) {
665 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
666 data.forcesOnlyOnEntitiesRow.end()) {
672 if (!data.forcesOnlyOnEntitiesCol.empty()) {
673 colIndices.resize(col_nb_dofs,
false);
675 col_indices = &colIndices[0];
677 VectorDofs::iterator dit = dofs.begin();
678 for (
int ii = 0; dit != dofs.end(); ++dit, ++ii) {
679 if (data.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
680 data.forcesOnlyOnEntitiesCol.end()) {
686 Mat
B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
687 : getFEMethod()->snes_B;
690 &*NN.data().begin(), ADD_VALUES);
702 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
707 if (col_type != MBVERTEX)
710 dataAtPts->faceRowData = &row_data;
711 CHKERR loopSideVolumes(sideFeName, *sideFe);
724 if (dataAtPts->faceRowData ==
nullptr)
727 if (row_type != MBVERTEX)
730 row_nb_dofs = dataAtPts->faceRowData->getIndices().size();
737 nb_gauss_pts = dataAtPts->faceRowData->getN().size1();
739 nb_base_fun_row = dataAtPts->faceRowData->getFieldData().size() / 3;
742 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col,
false);
746 CHKERR iNtegrate(*(dataAtPts->faceRowData), col_data);
749 CHKERR aSsemble(*(dataAtPts->faceRowData), col_data);
766 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2), &
m(
r + 1,
c + 0),
767 &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2), &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1),
771 auto t_w = getFTensor0IntegrationWeight();
773 auto t_normal = getFTensor1NormalsAtGaussPts();
775 auto t_inv_H = getFTensor2FromMat<3, 3>(dataAtPts->invHMat);
778 if (
auto arc_dof = dataAtPts->arcLengthDof.lock()) {
779 lambda = arc_dof->getFieldData();
782 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
784 double a = -0.5 * t_w * dAta.data.data.value1 *
lambda;
789 for (; bbc != nb_base_fun_col; ++bbc) {
794 for (; bbr != nb_base_fun_row; ++bbr) {
796 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
800 a * t_row_base * t_inv_H(
k,
i) * t_col_diff_base(
k) * t_normal(
j);
821 const int *row_indices = &*row_data.
getIndices().data().begin();
823 const int *col_indices = &*col_data.
getIndices().data().begin();
825 auto &data = *dataAtPts;
826 if (!data.forcesOnlyOnEntitiesRow.empty()) {
827 rowIndices.resize(row_nb_dofs,
false);
829 row_indices = &rowIndices[0];
831 VectorDofs::iterator dit = dofs.begin();
832 for (
int ii = 0; dit != dofs.end(); ++dit, ++ii) {
833 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
834 data.forcesOnlyOnEntitiesRow.end()) {
840 if (!data.forcesOnlyOnEntitiesCol.empty()) {
841 colIndices.resize(col_nb_dofs,
false);
843 col_indices = &colIndices[0];
845 VectorDofs::iterator dit = dofs.begin();
846 for (
int ii = 0; dit != dofs.end(); ++dit, ++ii) {
847 if (data.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
848 data.forcesOnlyOnEntitiesCol.end()) {
854 Mat
B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
855 : getFEMethod()->snes_B;
858 &*NN.data().begin(), ADD_VALUES);
877 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2), &
m(
r + 1,
c + 0),
878 &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2), &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1),
882 auto t_w = getFTensor0IntegrationWeight();
884 auto t_normal = getFTensor1NormalsAtGaussPts();
886 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hMat);
887 auto t_inv_H = getFTensor2FromMat<3, 3>(dataAtPts->invHMat);
892 if (
auto arc_dof = dataAtPts->arcLengthDof.lock()) {
893 lambda = arc_dof->getFieldData();
896 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
898 double a = -0.5 * t_w * dAta.data.data.value1 *
lambda;
903 for (; bbc != nb_base_fun_col; ++bbc) {
908 for (; bbr != nb_base_fun_row; ++bbr) {
910 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
914 t_assemble(
i,
j) -=
a * t_row_base * t_inv_H(
l,
j) *
915 t_col_diff_base(
m) * t_inv_H(
m,
i) * t_h(
k,
l) *
933 boost::ptr_vector<MethodForForceScaling> &methods_op,
bool ho_geometry)
936 F(
_F), dAta(data),
methodsOp(methods_op), hoGeometry(ho_geometry) {}
944 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
950 int nb_row_dofs = data.
getIndices().size() / rank;
955 EntityType fe_type = getNumeredEntFiniteElementPtr()->getEntType();
957 for (
unsigned int gg = 0; gg != data.
getN().size1(); ++gg) {
959 double val = getGaussPts()(2, gg);
961 if (hoGeometry || fe_type == MBQUAD) {
962 double area = cblas_dnrm2(3, &getNormalsAtGaussPts()(gg, 0), 1);
963 if (fe_type == MBTRI)
965 flux = dAta.data.data.value1 * area;
967 flux = dAta.data.data.value1 * getArea();
970 cblas_daxpy(nb_row_dofs, val *
flux, &data.
getN()(gg, 0), 1,
971 &*Nf.data().begin(), 1);
999 std::vector<double> mydata;
1002 for (
unsigned int ii = 0; ii != mydata.size(); ++ii) {
1003 force[ii] = mydata[ii];
1006 if (force.empty()) {
1011 const string name =
"Force";
1012 strncpy(
mapForce[ms_id].data.data.name, name.c_str(),
1013 name.size() > 5 ? 5 : name.size());
1015 std::sqrt(force[0] * force[0] + force[1] * force[1] + force[2] * force[2]);
1016 mapForce[ms_id].data.data.value1 = -magnitude;
1017 mapForce[ms_id].data.data.value2 = 0;
1019 force[0] / magnitude;
1021 force[1] / magnitude;
1023 force[2] / magnitude;
1024 mapForce[ms_id].data.data.value6 = 0;
1025 mapForce[ms_id].data.data.value7 = 0;
1026 mapForce[ms_id].data.data.value8 = 0;
1027 mapForce[ms_id].data.data.zero[0] = 0;
1028 mapForce[ms_id].data.data.zero[1] = 0;
1029 mapForce[ms_id].data.data.zero[2] = 0;
1030 mapForce[ms_id].data.data.zero2 = 0;
1060 &cubit_meshset_ptr);
1061 std::vector<double> mydata;
1064 for (
unsigned int ii = 0; ii != mydata.size(); ++ii) {
1065 pressure[ii] = mydata[ii];
1067 if (pressure.empty()) {
1070 const string name =
"Pressure";
1071 strncpy(
mapPressure[ms_id].data.data.name, name.c_str(),
1072 name.size() > 8 ? 8 : name.size());
1075 mapPressure[ms_id].data.data.value1 = pressure[0];
1093 const std::string x_field,
const std::string X_field,
1094 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1095 std::string side_fe_name,
Vec F, Mat aij,
int ms_id,
bool ho_geometry,
1104 &cubit_meshset_ptr);
1105 std::vector<double> mydata;
1108 for (
unsigned int ii = 0; ii != mydata.size(); ++ii) {
1109 pressure[ii] = mydata[ii];
1111 if (pressure.empty()) {
1114 const string name =
"Pressure";
1115 strncpy(
mapPressure[ms_id].data.data.name, name.c_str(),
1116 name.size() > 8 ? 8 : name.size());
1119 mapPressure[ms_id].data.data.value1 = pressure[0];
1136 x_field, X_field, data_at_pts, aij,
mapPressure[ms_id], ho_geometry));
1141 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideRhs =
1142 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(
mField);
1144 feMatSideRhs->getOpPtrVector().push_back(
1146 data_at_pts->getHMatPtr()));
1147 feMatSideRhs->getOpPtrVector().push_back(
1149 data_at_pts->getSmallhMatPtr()));
1151 feMatSideRhs->getOpPtrVector().push_back(
1155 X_field, data_at_pts, feMatSideRhs, side_fe_name, F,
mapPressure[ms_id],
1161 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideLhs_dx =
1162 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(
mField);
1164 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideLhs_dX =
1165 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(
mField);
1167 feMatSideLhs_dx->getOpPtrVector().push_back(
1169 data_at_pts->getHMatPtr()));
1170 feMatSideLhs_dx->getOpPtrVector().push_back(
1172 data_at_pts->getSmallhMatPtr()));
1174 feMatSideLhs_dx->getOpPtrVector().push_back(
1176 feMatSideLhs_dx->getOpPtrVector().push_back(
1178 X_field, x_field, data_at_pts, aij,
mapPressure[ms_id], ho_geometry));
1181 X_field, x_field, data_at_pts, feMatSideLhs_dx, side_fe_name, aij,
1184 feMatSideLhs_dX->getOpPtrVector().push_back(
1186 data_at_pts->getHMatPtr()));
1187 feMatSideLhs_dX->getOpPtrVector().push_back(
1189 data_at_pts->getSmallhMatPtr()));
1190 feMatSideLhs_dX->getOpPtrVector().push_back(
1192 feMatSideLhs_dX->getOpPtrVector().push_back(
1194 X_field, X_field, data_at_pts, aij,
mapPressure[ms_id], ho_geometry));
1198 X_field, X_field, data_at_pts, feMatSideLhs_dX, side_fe_name, aij,
1206 int ms_id,
bool ho_geometry) {
1213 std::vector<double> mydata;
1215 if (mydata.size() < 4)
1217 "Should be four block attributes but is %d", mydata.size());
1219 for (
unsigned int ii = 0; ii != 3; ++ii) {
1220 pressure_coeffs[ii] = mydata[ii];
1222 const double pressure_shift = mydata[3];
1227 boost::shared_ptr<MethodForAnalyticalForce> analytical_force_op(
1261 const std::string mesh_nodals_positions, Range *intersect_ptr) {
1272 mesh_nodals_positions);
1279 CHKERR m_field.
get_moab().get_entities_by_dimension(it->meshset, 2, tris,
1282 tris = intersect(tris, *intersect_ptr);
1293 mesh_nodals_positions);
1299 CHKERR m_field.
get_moab().get_entities_by_dimension(it->meshset, 2, tris,
1302 tris = intersect(tris, *intersect_ptr);
1308 const string block_set_force_name(
"FORCE");
1311 if (it->getName().compare(0, block_set_force_name.length(),
1312 block_set_force_name) == 0) {
1314 CHKERR m_field.
get_moab().get_entities_by_dimension(it->meshset, 2, tris,
1317 tris = intersect(tris, *intersect_ptr);
1323 const string block_set_pressure_name(
"PRESSURE");
1325 if (it->getName().compare(0, block_set_pressure_name.length(),
1326 block_set_pressure_name) == 0) {
1328 CHKERR m_field.
get_moab().get_entities_by_dimension(it->meshset, 2, tris,
1331 tris = intersect(tris, *intersect_ptr);
1338 const string block_set_linear_pressure_name(
"LINEAR_PRESSURE");
1340 if (it->getName().compare(0, block_set_linear_pressure_name.length(),
1341 block_set_linear_pressure_name) == 0) {
1343 CHKERR m_field.
get_moab().get_entities_by_dimension(it->meshset, 2, tris,
1346 tris = intersect(tris, *intersect_ptr);
1356 boost::ptr_map<std::string, NeumannForcesSurface> &neumann_forces,
Vec F,
1357 const std::string
field_name,
const std::string mesh_nodals_positions) {
1359 bool ho_geometry = m_field.
check_field(mesh_nodals_positions);
1362 std::string fe_name =
"FORCE_FE";
1364 auto &nf = neumann_forces.at(fe_name);
1365 auto &
fe = nf.getLoopFe();
1375 const string block_set_force_name(
"FORCE");
1377 if (it->getName().compare(0, block_set_force_name.length(),
1378 block_set_force_name) == 0) {
1387 std::string fe_name =
"PRESSURE_FE";
1390 auto &nf = neumann_forces.at(fe_name);
1391 auto &
fe = nf.getLoopFe();
1403 const string block_set_pressure_name(
"PRESSURE");
1405 if (it->getName().compare(0, block_set_pressure_name.length(),
1406 block_set_pressure_name) == 0) {
1413 const string block_set_linear_pressure_name(
"LINEAR_PRESSURE");
1415 if (it->getName().compare(0, block_set_linear_pressure_name.length(),
1416 block_set_linear_pressure_name) == 0) {
1428 const std::string mesh_nodals_positions) {
1437 mesh_nodals_positions);
1443 CHKERR m_field.
get_moab().get_entities_by_dimension(it->meshset, 2, tris,
1453 boost::ptr_map<std::string, NeumannForcesSurface> &neumann_forces,
Vec F,
1454 const std::string
field_name,
const std::string mesh_nodals_positions) {
1456 bool ho_geometry = m_field.
check_field(mesh_nodals_positions);
1459 fe_name =
"FLUX_FE";
1462 auto &nf = neumann_forces.at(fe_name);
1463 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.
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_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
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
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.
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.
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
double flux
impulse magnitude
constexpr auto field_name
FTensor::Index< 'm', 3 > m
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_vector< 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)
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)
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)
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.
NeumannForcesSurface(MoFEM::Interface &m_field)
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