9805 EntityHandle fe_ent = numeredEntFiniteElementPtr->getEnt();
9810 auto constrain_nodes = [&]() {
9812 std::vector<int> lambda_dofs(9, -1);
9813 std::vector<double> lambda_vals(9, 0);
9818 auto row_dofs = getRowDofsPtr();
9819 auto lo_uid_lambda = DofEntity::getLoFieldEntityUId(
9821 auto hi_uid_lambda = DofEntity::getHiFieldEntityUId(
9824 for (
auto it = row_dofs->lower_bound(lo_uid_lambda);
9825 it != row_dofs->upper_bound(hi_uid_lambda); ++it) {
9826 int side_number = it->get()->getSideNumberPtr()->side_number;
9827 if (side_number > 2)
9829 const int dof_idx = 3 * side_number + it->get()->getEntDofIdx();
9830 lambda_dofs[dof_idx] = it->get()->getPetscGlobalDofIdx();
9831 lambda_vals[dof_idx] = it->get()->getFieldData();
9836 "Only nine Lagrange multipliers can be on element (%d)", sum_dd);
9838 std::vector<int> positions_dofs(18, -1);
9839 std::vector<double> positions_vals(18, 0);
9840 std::vector<bool> positions_flags(18,
true);
9842 auto lo_uid_pos = DofEntity::getLoFieldEntityUId(
9844 auto hi_uid_pos = DofEntity::getHiFieldEntityUId(
9848 for (
auto it = row_dofs->lower_bound(lo_uid_pos);
9849 it != row_dofs->upper_bound(hi_uid_pos); ++it) {
9850 int side_number = it->get()->getSideNumberPtr()->side_number;
9851 int dof_idx = 3 * side_number + it->get()->getEntDofIdx();
9852 positions_dofs[dof_idx] = it->get()->getPetscGlobalDofIdx();
9853 positions_vals[dof_idx] = it->get()->getFieldData();
9854 Range vert =
Range(it->get()->getEnt(), it->get()->getEnt());
9856 positions_flags[dof_idx] =
false;
9860 const double beta = area *
aLpha;
9862 case CTX_SNESSETFUNCTION:
9863 CHKERR VecSetOption(snes_f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
9864 for (
int ii = 0; ii != 9; ++ii) {
9865 if (lambda_dofs[ii] == -1)
9867 double gap = positions_vals[0 + ii] - positions_vals[9 + ii];
9868 double val1 = beta * gap;
9869 CHKERR VecSetValue(snes_f, lambda_dofs[ii], val1, ADD_VALUES);
9870 double val2 = beta * lambda_vals[ii];
9871 if (positions_flags[0 + ii]) {
9872 CHKERR VecSetValue(snes_f, positions_dofs[0 + ii], +val2, ADD_VALUES);
9874 if (positions_flags[9 + ii]) {
9875 CHKERR VecSetValue(snes_f, positions_dofs[9 + ii], -val2, ADD_VALUES);
9879 case CTX_SNESSETJACOBIAN:
9880 for (
int ii = 0; ii != 9; ++ii) {
9881 if (lambda_dofs[ii] == -1)
9883 CHKERR MatSetValue(snes_B, lambda_dofs[ii], positions_dofs[0 + ii],
9884 +1 * beta, ADD_VALUES);
9885 CHKERR MatSetValue(snes_B, lambda_dofs[ii], positions_dofs[9 + ii],
9886 -1 * beta, ADD_VALUES);
9887 if (positions_flags[0 + ii]) {
9888 CHKERR MatSetValue(snes_B, positions_dofs[0 + ii], lambda_dofs[ii],
9889 +1 * beta, ADD_VALUES);
9891 if (positions_flags[9 + ii]) {
9892 CHKERR MatSetValue(snes_B, positions_dofs[9 + ii], lambda_dofs[ii],
9893 -1 * beta, ADD_VALUES);
9903 auto constrain_edges = [&]() {
9906 map<EntityHandle, EntityHandle> side_map;
9909 for (
auto s : {0, 1, 2}) {
9928 map<EntityHandle, std::vector<int>> side_map_lambda_dofs;
9929 map<EntityHandle, std::vector<double>> side_map_lambda_vals;
9930 map<EntityHandle, std::vector<int>> side_map_positions_dofs;
9931 map<EntityHandle, std::vector<double>> side_map_positions_vals;
9932 map<EntityHandle, double> side_map_positions_sign;
9934 for (
auto e : ents) {
9936 auto row_dofs = getRowDofsPtr();
9937 auto lo_uid_lambda = row_dofs->lower_bound(DofEntity::getLoFieldEntityUId(
9939 auto hi_uid_lambda = row_dofs->upper_bound(DofEntity::getHiFieldEntityUId(
9941 const int nb_lambda_dofs = std::distance(lo_uid_lambda, hi_uid_lambda);
9943 std::vector<int> lambda_dofs(nb_lambda_dofs, -1);
9944 std::vector<double> lambda_vals(nb_lambda_dofs, 0);
9946 for (
auto it = lo_uid_lambda; it != hi_uid_lambda; ++it) {
9947 const int dof_idx = it->get()->getEntDofIdx();
9948 lambda_dofs[dof_idx] = it->get()->getPetscGlobalDofIdx();
9949 lambda_vals[dof_idx] = it->get()->getFieldData();
9952 if (nb_lambda_dofs) {
9953 side_map_lambda_dofs[e] = lambda_dofs;
9954 side_map_lambda_vals[e] = lambda_vals;
9955 side_map_lambda_dofs[side_map[e]] = lambda_dofs;
9956 side_map_lambda_vals[side_map[e]] = lambda_vals;
9959 auto lo_uid_pos = row_dofs->lower_bound(DofEntity::getLoFieldEntityUId(
9961 auto hi_uid_pos = row_dofs->upper_bound(DofEntity::getHiFieldEntityUId(
9963 const int nb_dofs = std::distance(lo_uid_pos, hi_uid_pos);
9965 std::vector<int> positions_dofs(nb_dofs, -1);
9966 std::vector<double> positions_vals(nb_dofs, 0);
9967 side_map_positions_sign[e] = (!nb_lambda_dofs) ? 1 : -1;
9970 for (
auto it = lo_uid_pos; it != hi_uid_pos; ++it) {
9971 int dof_idx = it->get()->getEntDofIdx();
9972 positions_dofs[dof_idx] = it->get()->getPetscGlobalDofIdx();
9973 positions_vals[dof_idx] = it->get()->getFieldData();
9976 side_map_positions_dofs[e] = positions_dofs;
9977 side_map_positions_vals[e] = positions_vals;
9980 for (
auto m : side_map_positions_dofs) {
9983 if (side_map_lambda_dofs.find(e) != side_map_lambda_dofs.end()) {
9985 auto &lambda_dofs = side_map_lambda_dofs.at(e);
9986 auto &lambda_vals = side_map_lambda_vals.at(e);
9987 auto &positions_dofs = side_map_positions_dofs.at(e);
9988 auto &positions_vals = side_map_positions_vals.at(e);
9989 auto sign = side_map_positions_sign.at(e);
9991 const double beta = area *
aLpha;
9994 case CTX_SNESSETFUNCTION:
9995 for (
int ii = 0; ii != lambda_dofs.size(); ++ii) {
9996 double val1 = beta * positions_vals[ii] *
sign;
9997 CHKERR VecSetValue(snes_f, lambda_dofs[ii], val1, ADD_VALUES);
9999 for (
int ii = 0; ii != positions_dofs.size(); ++ii) {
10000 double val2 =
sign * beta * lambda_vals[ii];
10001 CHKERR VecSetValue(snes_f, positions_dofs[ii], +val2, ADD_VALUES);
10004 case CTX_SNESSETJACOBIAN:
10005 for (
int ii = 0; ii != lambda_dofs.size(); ++ii) {
10006 CHKERR MatSetValue(snes_B, lambda_dofs[ii], positions_dofs[ii],
10007 beta *
sign, ADD_VALUES);
10009 for (
int ii = 0; ii != positions_dofs.size(); ++ii) {
10010 CHKERR MatSetValue(snes_B, positions_dofs[ii], lambda_dofs[ii],
10011 beta *
sign, ADD_VALUES);
10023 CHKERR constrain_nodes();
10024 CHKERR constrain_edges();