ElasticTie.hpp

Go to the documentation of this file.

/**
 * @file ElasticTie.hpp
 *
 * @brief TIE constraint support for the elastic tutorial
 *
 * @copyright Copyright (c) 2025
 *
 */
 
#ifndef __ELASTIC_TIE_HPP__
#define __ELASTIC_TIE_HPP__
 
#include <regex>
#include <quad.h>
 
namespace ElasticTie {
 
FTENSOR_INDEX(SPACE_DIM, i);
 
// we know this because we created the prisms in insert_constraint_prisms.cpp
constexpr int SLAVE_FACE_SIDE = 3;
constexpr int MASTER_FACE_SIDE = 4;
 
struct CommonData {
  bool hasTieConstraints = false;
  std::string tiePrismFeName = "TIE_PRISM_FE";
  std::string dispFieldName = "U";
  FieldApproximationBase dispFieldBase = AINSWORTH_LEGENDRE_BASE;
  int dispFieldOrder = 1;
  int tieRefLevel = 1;
  boost::shared_ptr<Range> tieMasterFaces = nullptr;
  boost::shared_ptr<Range> tieSlaveFaces = nullptr;
  boost::shared_ptr<Range> tiePrisms = nullptr;
  boost::shared_ptr<MatrixDouble> dispAtGaussPts =
      boost::make_shared<MatrixDouble>();
  boost::shared_ptr<MatrixDouble> lambdaAtGaussPts =
      boost::make_shared<MatrixDouble>();
  boost::shared_ptr<FlatPrismElementForcesAndSourcesCore> tiePrismRhsFE =
      nullptr;
  boost::shared_ptr<FlatPrismElementForcesAndSourcesCore> tiePrismLhsFE =
      nullptr;
};
 
 
inline bool isInsideReferenceTriangle(const double xi, const double eta,
                                      const double tol = 1e-8) {
  return xi >= -tol && eta >= -tol && xi + eta <= 1 + tol;
}
 
inline bool isOnPrismTriFace(const int side, const EntityType type,
                             const int prism_face_side) {
  switch (type) {
  case MBVERTEX:
    return prism_face_side == SLAVE_FACE_SIDE ? (side >= 0 && side < 3)
                                              : (side >= 3 && side < 6);
  case MBEDGE:
    return prism_face_side == SLAVE_FACE_SIDE ? (side >= 0 && side < 3)
                                              : (side >= 6 && side < 9);
  case MBTRI:
    return side == prism_face_side;
  default:
    return false;
  }
}
 
inline bool isOnPrismTieTrace(const int side, const EntityType type) {
  return isOnPrismTriFace(side, type, SLAVE_FACE_SIDE) ||
         isOnPrismTriFace(side, type, MASTER_FACE_SIDE);
}
 
inline bool isOnSlaveTieRowBase(const int rr, const int row_side,
                                const EntityType row_type) {
  switch (row_type) {
  case MBVERTEX:
    return rr < 3;
  case MBEDGE:
    return row_side >= 0 && row_side < 3;
  case MBTRI:
    return row_side == SLAVE_FACE_SIDE;
  default:
    return false;
  }
}
 
struct SetIntegrationOnTiePrismFaces {
 
  using FunRule = boost::function<int(int)>;
  FunRule funRule = [](int p) { return 2 * p + 1; };
  int tieRefLevel = 1;
 
  SetIntegrationOnTiePrismFaces() = delete;
 
  SetIntegrationOnTiePrismFaces(int tie_ref_level = 1)
      : tieRefLevel(tie_ref_level) {}
 
  SetIntegrationOnTiePrismFaces(FunRule fun_rule, int tie_ref_level = 1)
      : funRule(fun_rule), tieRefLevel(tie_ref_level) {}
 
  MoFEMErrorCode operator()(ForcesAndSourcesCore *fe_raw_ptr, int order_row,
                            int order_col, int order_data) {
    MoFEMFunctionBegin;
 
    auto fe_ptr = static_cast<FlatPrismElementForcesAndSourcesCore *>(fe_raw_ptr);
    const int rule = funRule(std::max(order_data, std::max(order_row, order_col)));
 
    if (tieRefLevel < 0) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "-tie_ref_level has to be >= 0");
    }
 
    if (rule < 0 || rule >= QUAD_2D_TABLE_SIZE) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Triangle quadrature rule %d is not available for TIE prism face "
              "integration",
              rule);
    }
 
    if (QUAD_2D_TABLE[rule]->dim != 2) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Expected 2D quadrature for TIE prism face integration");
    }
 
    if (QUAD_2D_TABLE[rule]->order < rule) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Wrong quadrature order %d < %d for TIE prism face integration",
              QUAD_2D_TABLE[rule]->order, rule);
    }
 
    CHKERR setBaseQuadrature(*fe_ptr, rule);
 
    if (tieRefLevel > 0) {
      const auto cache_key = std::make_pair(rule, tieRefLevel);
      auto it = mapRefCoords.find(cache_key);
      if (it == mapRefCoords.end()) {
        MatrixDouble ref_gauss_pts;
        CHKERR refineQuadrature(*fe_ptr, tieRefLevel, ref_gauss_pts);
        mapRefCoords[cache_key].swap(ref_gauss_pts);
        it = mapRefCoords.find(cache_key);
      }
      fe_ptr->gaussPts = it->second;
    }
 
    const MatrixDouble slave_gauss_pts = fe_ptr->gaussPts;
    const auto nb_gauss_pts = slave_gauss_pts.size2();
 
    auto &moab = fe_ptr->mField.get_moab();
    auto get_prism_face_coords = [&](const int side, std::array<double, 9> &coords) {
      MoFEMFunctionBegin;
      const EntityHandle prism = fe_ptr->getFEEntityHandle();
      const EntityHandle *prism_conn = nullptr;
      int num_prism_nodes = 0;
      CHKERR moab.get_connectivity(prism, prism_conn, num_prism_nodes, true);
      if (num_prism_nodes != 6) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "TIE prism is expected to have 6 nodes");
      }
 
      std::array<EntityHandle, 3> face_conn;
      switch (side) {
      case SLAVE_FACE_SIDE:
        face_conn = {prism_conn[0], prism_conn[1], prism_conn[2]};
        break;
      case MASTER_FACE_SIDE:
        face_conn = {prism_conn[3], prism_conn[4], prism_conn[5]};
        break;
      default:
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Unsupported TIE prism triangle side %d", side);
      }
 
      CHKERR moab.get_coords(face_conn.data(), 3, coords.data());
      MoFEMFunctionReturn(0);
    };
 
    std::array<double, 9> slave_coords; // nodal coords
    std::array<double, 9> master_coords; // nodal coords
    CHKERR get_prism_face_coords(SLAVE_FACE_SIDE, slave_coords);
    CHKERR get_prism_face_coords(MASTER_FACE_SIDE, master_coords);
 
    MatrixDouble slave_global_coords(nb_gauss_pts, 3, false); // gp coords
    MatrixDouble master_global_coords(nb_gauss_pts, 3, false); //gp coords
    MatrixDouble slave_shape(nb_gauss_pts, 3, false);
    CHKERR ShapeMBTRI(&slave_shape(0, 0), &slave_gauss_pts(0, 0),
                      &slave_gauss_pts(1, 0), nb_gauss_pts);
    for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
      for (int dd = 0; dd != 3; ++dd) {
        slave_global_coords(gg, dd) =
            slave_shape(gg, 0) * slave_coords[0 + dd] +
            slave_shape(gg, 1) * slave_coords[3 + dd] +
            slave_shape(gg, 2) * slave_coords[6 + dd];
      }
    }
 
    CHKERR projectSlavePointsToMasterPlane(master_coords.data(),
                                           slave_global_coords,
                                           master_global_coords);
 
    MatrixDouble master_local_coords(nb_gauss_pts, 2, false);
    CHKERR Tools::getLocalCoordinatesOnReferenceThreeNodeTri(
        master_coords.data(), &master_global_coords(0, 0),
        master_global_coords.size1(), &master_local_coords(0, 0));
 
    MatrixDouble combined_gauss_pts(3, 2 * nb_gauss_pts, false);
    for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
      const bool valid_master_point =
          isInsideReferenceTriangle(master_local_coords(gg, 0),
                                    master_local_coords(gg, 1));
      const double paired_weight = valid_master_point ? slave_gauss_pts(2, gg) : 0.;
 
      combined_gauss_pts(0, gg) = slave_gauss_pts(0, gg);
      combined_gauss_pts(1, gg) = slave_gauss_pts(1, gg);
      combined_gauss_pts(2, gg) = paired_weight;
 
      combined_gauss_pts(0, gg + nb_gauss_pts) = master_local_coords(gg, 0);
      combined_gauss_pts(1, gg + nb_gauss_pts) = master_local_coords(gg, 1);
      combined_gauss_pts(2, gg + nb_gauss_pts) = paired_weight;
    }
    fe_ptr->gaussPts.swap(combined_gauss_pts);
 
    MoFEMFunctionReturn(0);
  }
 
private:
  static MoFEMErrorCode setBaseQuadrature(
      FlatPrismElementForcesAndSourcesCore &fe, const int rule) {
    MoFEMFunctionBegin;
    const size_t nb_gauss_pts = QUAD_2D_TABLE[rule]->npoints;
    fe.gaussPts.resize(3, nb_gauss_pts, false);
    cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[1], 3,
                &fe.gaussPts(0, 0), 1);
    cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[2], 3,
                &fe.gaussPts(1, 0), 1);
    cblas_dcopy(nb_gauss_pts, QUAD_2D_TABLE[rule]->weights, 1,
                &fe.gaussPts(2, 0), 1);
    MoFEMFunctionReturn(0);
  }
 
  static MoFEMErrorCode projectPointToMasterPlane(const double *master_coords,
                                                  const double *point,
                                                  double *projected_point) {
    MoFEMFunctionBegin;
 
    // For the infinite master plane, the closest point is the orthogonal
    // projection along the plane normal.
    const double ax = master_coords[0];
    const double ay = master_coords[1];
    const double az = master_coords[2];
    const double bx = master_coords[3];
    const double by = master_coords[4];
    const double bz = master_coords[5];
    const double cx = master_coords[6];
    const double cy = master_coords[7];
    const double cz = master_coords[8];
 
    const double abx = bx - ax;
    const double aby = by - ay;
    const double abz = bz - az;
    const double acx = cx - ax;
    const double acy = cy - ay;
    const double acz = cz - az;
 
    const double nx = aby * acz - abz * acy;
    const double ny = abz * acx - abx * acz;
    const double nz = abx * acy - aby * acx;
    const double n_norm_sq = nx * nx + ny * ny + nz * nz;
    if (n_norm_sq <= std::numeric_limits<double>::epsilon()) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Degenerated master triangle plane");
    }
 
    const double px = point[0] - ax;
    const double py = point[1] - ay;
    const double pz = point[2] - az;
    const double signed_distance = (px * nx + py * ny + pz * nz) / n_norm_sq;
    projected_point[0] = point[0] - signed_distance * nx;
    projected_point[1] = point[1] - signed_distance * ny;
    projected_point[2] = point[2] - signed_distance * nz;
 
    MoFEMFunctionReturn(0);
  }
 
  static MoFEMErrorCode projectSlavePointsToMasterPlane(
      const double *master_coords, const MatrixDouble &slave_global_coords,
      MatrixDouble &master_global_coords) {
    MoFEMFunctionBegin;
 
    // The weak form is integrated on the slave surface. Each slave Gauss point
    // is paired with its orthogonal projection on the infinite plane defined
    // by the master triangle.
    master_global_coords.resize(slave_global_coords.size1(), 3, false);
    for (size_t gg = 0; gg != slave_global_coords.size1(); ++gg) {
      CHKERR projectPointToMasterPlane(master_coords,
                                       &slave_global_coords(gg, 0),
                                       &master_global_coords(gg, 0));
    }
 
    MoFEMFunctionReturn(0);
  }
 
  static MoFEMErrorCode refineQuadrature(
      FlatPrismElementForcesAndSourcesCore &fe, const int refinement_levels,
      MatrixDouble &ref_gauss_pts) {
    MoFEMFunctionBegin;
 
    constexpr int num_nodes = 3;
    moab::Core moab_ref;
    double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0};
    EntityHandle nodes[num_nodes];
    for (int nn = 0; nn != num_nodes; ++nn)
      CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
 
    EntityHandle tri;
    CHKERR moab_ref.create_element(MBTRI, nodes, num_nodes, tri);
    MoFEM::CoreTmp<-1> mofem_ref_core(moab_ref, PETSC_COMM_SELF, -2);
    MoFEM::Interface &m_field_ref = mofem_ref_core;
 
    {
      Range tris(tri, tri);
      Range edges;
      // Force MOAB to materialise the triangle edges before the bit-ref
      // database is seeded. The working reference-refinement snippets all do
      // this, and the side-number lookup below verifies that MOAB sees a
      // consistent local edge numbering on the temporary reference triangle.
      CHKERR m_field_ref.get_moab().get_adjacencies(
          tris, 1, true, edges, moab::Interface::UNION);
      for (auto edge : edges) {
        int side_number = -1, sense = 0, offset = 0;
        CHKERR moab_ref.side_number(tri, edge, side_number, sense, offset);
      }
      CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevel(
          tris, BitRefLevel().set(0), false, VERBOSE);
    }
 
    auto *m_ref = m_field_ref.getInterface<MeshRefinement>();
    for (int ll = 0; ll != refinement_levels; ++ll) {
      Range tris;
      CHKERR m_field_ref.getInterface<BitRefManager>()
          ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
                                         BitRefLevel().set(), MBTRI, tris);
 
      Range edges;
      CHKERR moab_ref.get_adjacencies(tris, 1, true, edges,
                                      moab::Interface::UNION);
 
      const BitRefLevel child_bit = BitRefLevel().set(ll + 1);
      CHKERR m_ref->addVerticesInTheMiddleOfEdges(edges, child_bit);
      CHKERR m_ref->refineTris(tris, child_bit, QUIET, false);
    }
 
    Range tris;
    CHKERR m_field_ref.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
        BitRefLevel().set(refinement_levels), BitRefLevel().set(), MBTRI, tris);
 
    MatrixDouble ref_coords(tris.size(), 9, false);
    int tt = 0;
    for (auto tit = tris.begin(); tit != tris.end(); ++tit, ++tt) {
      int num_conn;
      const EntityHandle *conn;
      CHKERR moab_ref.get_connectivity(*tit, conn, num_conn, false);
      CHKERR moab_ref.get_coords(conn, num_conn, &ref_coords(tt, 0));
    }
 
    const size_t nb_gauss_pts = fe.gaussPts.size2();
    MatrixDouble shape_n(nb_gauss_pts, 3, false);
    CHKERR ShapeMBTRI(&shape_n(0, 0), &fe.gaussPts(0, 0), &fe.gaussPts(1, 0),
                      nb_gauss_pts);
    ref_gauss_pts.resize(3, nb_gauss_pts * ref_coords.size1(), false);
    int gg = 0;
    for (size_t rr = 0; rr != ref_coords.size1(); ++rr) {
      double *tri_coords = &ref_coords(rr, 0);
      FTensor::Tensor1<double, 3> t_normal;
      CHKERR Tools::getTriNormal(tri_coords, &t_normal(0));
      const double det = t_normal.l2();
      for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
        for (int dd = 0; dd != 2; ++dd) {
          ref_gauss_pts(dd, gg) =
              shape_n(ggg, 0) * tri_coords[3 * 0 + dd] +
              shape_n(ggg, 1) * tri_coords[3 * 1 + dd] +
              shape_n(ggg, 2) * tri_coords[3 * 2 + dd];
        }
        ref_gauss_pts(2, gg) = fe.gaussPts(2, ggg) * det;
      }
    }
 
    MoFEMFunctionReturn(0);
  }
 
  static inline std::map<std::pair<int, int>, MatrixDouble> mapRefCoords;
};
 
template <int Tensor_Dim>
struct OpCalculateTraceVectorFieldValues
    : public OpCalculateVectorFieldValues_General<Tensor_Dim, MatrixDouble> {
 
  using Base = OpCalculateVectorFieldValues_General<Tensor_Dim, MatrixDouble>;
  using Base::Base;
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data) override {
    MoFEMFunctionBeginHot;
    if (!isOnPrismTieTrace(side, type))
      MoFEMFunctionReturnHot(0);
    CHKERR Base::doWork(side, type, data);
    MoFEMFunctionReturnHot(0);
  }
};
 
using TiePrismRhsBase = FormsIntegrators<
    FlatPrismElementForcesAndSourcesCore::UserDataOperator>::Assembly<A>::OpBase;
 
struct OpTiePrismRhs : public TiePrismRhsBase {
 
  OpTiePrismRhs(boost::shared_ptr<MatrixDouble> disp_ptr)
      : TiePrismRhsBase("LAMBDA", "LAMBDA", OPROW, nullptr),
        dispPtr(disp_ptr) {
    if (!dispPtr)
      THROW_MESSAGE("Pointers are not set");
  }
 
  inline double getSlaveFaceArea() { return getAreaF3(); }
 
protected:
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data) override {
    MoFEMFunctionBeginHot;
 
    const auto *fe =
        dynamic_cast<const FlatPrismElementForcesAndSourcesCore *>(ptrFE);
    if (!fe)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "TIE prism RHS operator is attached to a wrong finite element");
 
    const auto nb_gauss_pts = getGaussPts().size2();
    if (!nb_gauss_pts || nb_gauss_pts % 2) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Expected concatenated [slave|master] TIE Gauss points");
    }
    const auto nb_slave_gauss = nb_gauss_pts / 2;
    if (dispPtr->size1() != nb_gauss_pts) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "TIE displacement values at Gauss points were not prepared "
              "before RHS assembly");
    }
 
    if (nbRows % SPACE_DIM)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Unexpected number of LAMBDA dofs on TIE prism rows");
 
    auto &row_n = row_data.getN();
    if (!row_n.size1() || !row_n.size2())
      MoFEMFunctionReturnHot(0);
    if (row_n.size1() != nb_gauss_pts)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Unexpected number of Gauss points on TIE multiplier block");
 
    const int nb_row_base = nbRows / SPACE_DIM;
    if (static_cast<size_t>(nb_row_base) > row_n.size2())
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "More TIE multiplier bases than shape functions");
    if (dispPtr->size2() != SPACE_DIM)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Unexpected TIE displacement matrix size: size2=%d expected=%d",
              static_cast<int>(dispPtr->size2()), SPACE_DIM);
    if (locF.size() != static_cast<size_t>(nbRows))
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Unexpected TIE RHS block size: rowSide=%d rowType=%d nbRows=%d "
              "locF.size()=%d row_n=(%d,%d) disp=(%d,%d)",
              rowSide, rowType, nbRows, static_cast<int>(locF.size()),
              static_cast<int>(row_n.size1()), static_cast<int>(row_n.size2()),
              static_cast<int>(dispPtr->size1()),
              static_cast<int>(dispPtr->size2()));
 
    auto t_w = getFTensor0IntegrationWeight();
    for (size_t gg = 0; gg != nb_slave_gauss; ++gg) {
      const double slave_alpha = getSlaveFaceArea() * t_w;
      auto t_row_base = row_data.getFTensor0N(gg, 0);
      for (int rr = 0; rr != nb_row_base; ++rr) {
        const double row_base = t_row_base;
        if (row_base != 0.) {
          auto t_row_rhs = getFTensor1FromPtr<SPACE_DIM>(&locF[SPACE_DIM * rr]);
          for (int dd = 0; dd != SPACE_DIM; ++dd) {
            t_row_rhs(dd) += row_base * slave_alpha *
                             ((*dispPtr)(gg, dd) -
                              (*dispPtr)(gg + nb_slave_gauss, dd));
          }
        }
        ++t_row_base;
      }
 
      ++t_w;
    }
 
    MoFEMFunctionReturnHot(0);
  }
 
private:
  boost::shared_ptr<MatrixDouble> dispPtr;
};
 
struct OpTiePrismRhsU : public TiePrismRhsBase {
 
  OpTiePrismRhsU(const std::string &field_name,
                 boost::shared_ptr<MatrixDouble> lambda_ptr)
      : TiePrismRhsBase(field_name, field_name, OPROW, nullptr),
        lambdaPtr(lambda_ptr) {
    if (!lambdaPtr)
      THROW_MESSAGE("Pointers are not set");
  }
 
  inline double getSlaveFaceArea() { return getAreaF3(); }
 
protected:
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data) override {
    MoFEMFunctionBeginHot;
 
    const auto *fe =
        dynamic_cast<const FlatPrismElementForcesAndSourcesCore *>(ptrFE);
    if (!fe)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "TIE prism RHS operator is attached to a wrong finite element");
 
    const auto nb_gauss_pts = getGaussPts().size2();
    if (!nb_gauss_pts || nb_gauss_pts % 2) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Expected concatenated [slave|master] TIE Gauss points");
    }
    const auto nb_slave_gauss = nb_gauss_pts / 2;
    if (lambdaPtr->size1() != nb_gauss_pts) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "TIE multiplier values at Gauss points were not prepared "
              "before RHS assembly");
    }
 
    if (nbRows % SPACE_DIM)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Unexpected number of TIE prism vector dofs on rows");
 
    if (!isOnPrismTieTrace(rowSide, rowType))
      MoFEMFunctionReturnHot(0);
 
    auto &row_n = row_data.getN();
    if (!row_n.size1() || !row_n.size2())
      MoFEMFunctionReturnHot(0);
    if (row_n.size1() != nb_gauss_pts)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Unexpected number of Gauss points on TIE displacement block");
 
    const int nb_row_base = nbRows / SPACE_DIM;
    if (static_cast<size_t>(nb_row_base) > row_n.size2())
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "More TIE displacement bases than shape functions");
    if (lambdaPtr->size2() != SPACE_DIM)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Unexpected TIE multiplier matrix size: size2=%d expected=%d",
              static_cast<int>(lambdaPtr->size2()), SPACE_DIM);
    if (locF.size() != static_cast<size_t>(nbRows))
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Unexpected TIE U RHS block size: rowSide=%d rowType=%d "
              "nbRows=%d locF.size()=%d row_n=(%d,%d) lambda=(%d,%d)",
              rowSide, rowType, nbRows, static_cast<int>(locF.size()),
              static_cast<int>(row_n.size1()), static_cast<int>(row_n.size2()),
              static_cast<int>(lambdaPtr->size1()),
              static_cast<int>(lambdaPtr->size2()));
 
    auto t_w = getFTensor0IntegrationWeight();
    for (size_t gg = 0; gg != nb_slave_gauss; ++gg) {
      const double slave_alpha = getSlaveFaceArea() * t_w;
      for (int rr = 0; rr != nb_row_base; ++rr) {
        const bool base_on_slave = isOnSlaveTieRowBase(rr, rowSide, rowType);
        const double row_sign = base_on_slave ? 1.0 : -1.0;
        const auto row_gg = base_on_slave ? gg : gg + nb_slave_gauss;
        if (row_gg >= row_n.size1())
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "TIE U RHS Gauss index out of range: rowSide=%d rowType=%d "
                  "rr=%d row_gg=%d row_n.size1()=%d nb_slave_gauss=%d",
                  rowSide, rowType, rr, static_cast<int>(row_gg),
                  static_cast<int>(row_n.size1()),
                  static_cast<int>(nb_slave_gauss));
        auto t_row_base = row_data.getFTensor0N(row_gg, rr);
        const double row_base = t_row_base;
        if (row_base != 0.) {
          auto t_row_rhs = getFTensor1FromPtr<SPACE_DIM>(&locF[SPACE_DIM * rr]);
          for (int dd = 0; dd != SPACE_DIM; ++dd)
            t_row_rhs(dd) +=
                row_sign * row_base * slave_alpha * (*lambdaPtr)(gg, dd);
        }
        ++t_row_base;
      }
 
      ++t_w;
    }
 
    MoFEMFunctionReturnHot(0);
  }
 
private:
  boost::shared_ptr<MatrixDouble> lambdaPtr;
};
 
using TiePrismLhsBase = FormsIntegrators<
    FlatPrismElementForcesAndSourcesCore::UserDataOperator>::Assembly<A>::OpBase;
 
struct OpTiePrismLhs : public TiePrismLhsBase {
 
  OpTiePrismLhs(const std::string &row_field_name,
                const std::string &col_field_name,
                const bool assemble_transpose = false)
      : TiePrismLhsBase(row_field_name, col_field_name, OPROWCOL, nullptr),
        rowFieldName(row_field_name), colFieldName(col_field_name) {
    sYmm = false;
    this->assembleTranspose = assemble_transpose;
  }
 
  inline double getSlaveFaceArea() { return getAreaF3(); }
 
protected:
  std::string rowFieldName;
  std::string colFieldName;
 
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
                           EntitiesFieldData::EntData &col_data) override {
    MoFEMFunctionBeginHot;
 
    const auto *fe =
        dynamic_cast<const FlatPrismElementForcesAndSourcesCore *>(ptrFE);
    if (!fe)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "TIE prism operator is attached to a wrong finite element");
 
    if (nbRows % SPACE_DIM || nbCols % SPACE_DIM)
      SETERRQ(
          PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
          "Unexpected number of TIE prism vector dofs on rows or columns: "
          "rowField=%s colField=%s rowSide=%d colSide=%d rowType=%d colType=%d "
          "nbRows=%d nbCols=%d",
          rowFieldName.c_str(), colFieldName.c_str(), rowSide, colSide,
          rowType, colType, nbRows, nbCols);
 
    if (colFieldName != "LAMBDA")
      MoFEMFunctionReturnHot(0);
 
    const auto nb_gauss_pts = getGaussPts().size2();
    if (!nb_gauss_pts || nb_gauss_pts % 2) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Expected concatenated [slave|master] TIE Gauss points");
    }
    const auto nb_slave_gauss = nb_gauss_pts / 2;
    const int nb_row_base = nbRows / SPACE_DIM;
    const int nb_col_base = nbCols / SPACE_DIM;
 
    auto &col_n = col_data.getN();
    if (!col_n.size1() || !col_n.size2())
      MoFEMFunctionReturnHot(0);
 
    auto &row_n = row_data.getN();
    if (!row_n.size1() || !row_n.size2())
      MoFEMFunctionReturnHot(0);
 
    if (row_n.size1() != nb_gauss_pts || col_n.size1() != nb_gauss_pts)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Unexpected number of Gauss points on TIE prism blocks");
 
    if (static_cast<size_t>(nb_row_base) > row_n.size2() ||
        static_cast<size_t>(nb_col_base) > col_n.size2())
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "More TIE prism bases than shape functions");
    if (locMat.size1() != static_cast<size_t>(nbRows) ||
        locMat.size2() != static_cast<size_t>(nbCols))
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Unexpected TIE LHS block size: rowField=%s colField=%s "
              "rowSide=%d colSide=%d rowType=%d colType=%d nbRows=%d "
              "nbCols=%d locMat=(%d,%d) row_n=(%d,%d) col_n=(%d,%d)",
              rowFieldName.c_str(), colFieldName.c_str(), rowSide, colSide,
              rowType, colType, nbRows, nbCols,
              static_cast<int>(locMat.size1()), static_cast<int>(locMat.size2()),
              static_cast<int>(row_n.size1()), static_cast<int>(row_n.size2()),
              static_cast<int>(col_n.size1()), static_cast<int>(col_n.size2()));
 
    if (!isOnPrismTieTrace(rowSide, rowType))
      MoFEMFunctionReturnHot(0);
 
    auto t_w = getFTensor0IntegrationWeight();
    for (size_t gg = 0; gg != nb_slave_gauss; ++gg) {
      const double slave_alpha = getSlaveFaceArea() * t_w;
      for (int rr = 0; rr != nb_row_base; ++rr) {
        const bool base_on_slave = isOnSlaveTieRowBase(rr, rowSide, rowType);
        const double row_sign = base_on_slave ? 1.0 : -1.0;
        const auto row_gg = base_on_slave ? gg : gg + nb_slave_gauss;
        if (row_gg >= row_n.size1())
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "TIE LHS Gauss index out of range: rowField=%s colField=%s "
                  "rowSide=%d rowType=%d rr=%d row_gg=%d row_n.size1()=%d "
                  "nb_slave_gauss=%d",
                  rowFieldName.c_str(), colFieldName.c_str(), rowSide, rowType,
                  rr, static_cast<int>(row_gg), static_cast<int>(row_n.size1()),
                  static_cast<int>(nb_slave_gauss));
        auto t_row_base = row_data.getFTensor0N(row_gg, rr);
        auto t_col_base = col_data.getFTensor0N(gg, 0);
        for (int cc = 0; cc != nb_col_base; ++cc) {
          const double value = row_sign * slave_alpha * t_row_base * t_col_base;
          for (int dd = 0; dd != SPACE_DIM; ++dd)
            locMat(SPACE_DIM * rr + dd, SPACE_DIM * cc + dd) += value;
          ++t_col_base;
        }
        ++t_row_base;
      }
      ++t_w;
    }
 
    MoFEMFunctionReturnHot(0);
  }
};
 
inline MoFEMErrorCode logTieDisplacementNorms(
    MoFEM::Interface &mField, const boost::shared_ptr<CommonData> &tie_data,
    const std::string &field_name) {
  MoFEMFunctionBegin;
 
  if (!tie_data || !tie_data->hasTieConstraints || !tie_data->tieSlaveFaces ||
      !tie_data->tieMasterFaces) {
    return 0;
  }
 
  auto &moab = mField.get_moab();
  auto comm_interface = mField.getInterface<CommInterface>();
  auto *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
 
  auto closure_with_faces = [&](const Range &faces, Range &closure) {
    MoFEMFunctionBegin;
    closure = faces;
    Range edges, verts;
    CHKERR moab.get_adjacencies(faces, 1, true, edges, moab::Interface::UNION);
    CHKERR moab.get_adjacencies(faces, 0, true, verts, moab::Interface::UNION);
    closure.merge(edges);
    closure.merge(verts);
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_norms = [&](const Range &ents, double &l2_norm, double &linf,
                             unsigned int &nb_dofs) {
    MoFEMFunctionBegin;
    double sum2 = 0;
    double max_abs = 0;
    unsigned int count = 0;
    for (auto ent : ents) {
      if (pcomm) {
        unsigned char pstatus;
        CHKERR moab.tag_get_data(pcomm->pstatus_tag(), &ent, 1, &pstatus);
        if (pstatus & PSTATUS_NOT_OWNED)
          continue;
      }
      for (_IT_GET_DOFS_FIELD_BY_NAME_AND_ENT_FOR_LOOP_(mField, field_name, ent,
                                                        dit)) {
        const double value = (*dit)->getFieldData();
        sum2 += value * value;
        max_abs = std::max(max_abs, std::abs(value));
        ++count;
      }
    }
    double global_sum2 = 0;
    double global_max_abs = 0;
    unsigned long long global_count = 0;
    const auto local_count = static_cast<unsigned long long>(count);
    MPI_Allreduce(&sum2, &global_sum2, 1, MPI_DOUBLE, MPI_SUM,
                  mField.get_comm());
    MPI_Allreduce(&max_abs, &global_max_abs, 1, MPI_DOUBLE, MPI_MAX,
                  mField.get_comm());
    MPI_Allreduce(&local_count, &global_count, 1, MPI_UNSIGNED_LONG_LONG,
                  MPI_SUM, mField.get_comm());
    l2_norm = std::sqrt(global_sum2);
    linf = global_max_abs;
    nb_dofs = static_cast<unsigned int>(global_count);
    MoFEMFunctionReturn(0);
  };
 
  Range slave_closure, master_closure;
  CHKERR closure_with_faces(*tie_data->tieSlaveFaces, slave_closure);
  CHKERR closure_with_faces(*tie_data->tieMasterFaces, master_closure);
  CHKERR comm_interface->synchroniseEntities(slave_closure);
  CHKERR comm_interface->synchroniseEntities(master_closure);
 
  double slave_l2 = 0, slave_linf = 0, master_l2 = 0, master_linf = 0;
  unsigned int slave_nb_dofs = 0, master_nb_dofs = 0;
  CHKERR calculate_norms(slave_closure, slave_l2, slave_linf, slave_nb_dofs);
  CHKERR calculate_norms(master_closure, master_l2, master_linf, master_nb_dofs);
 
  MOFEM_LOG("WORLD", Sev::inform)
      << "TIE solved displacement norms on slave face closure: dofs="
      << slave_nb_dofs << " l2=" << slave_l2 << " linf=" << slave_linf;
  MOFEM_LOG("WORLD", Sev::inform)
      << "TIE solved displacement norms on master face closure: dofs="
      << master_nb_dofs << " l2=" << master_l2 << " linf=" << master_linf;
 
  MoFEMFunctionReturn(0);
}
 
inline MoFEMErrorCode initializeTieConstraints(
    MoFEM::Interface &mField, boost::shared_ptr<CommonData> &tie_data,
    const std::string &field_name, FieldApproximationBase base, int order) {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto meshset_mng = mField.getInterface<MeshsetsManager>();
  auto comm_interface = mField.getInterface<CommInterface>();
 
  auto tie_master =
      meshset_mng->getCubitMeshsetPtr(std::regex("^TIE_MASTER.*"));
  auto tie_slave = meshset_mng->getCubitMeshsetPtr(std::regex("^TIE_SLAVE.*"));
 
  const bool have_master_meshset = !tie_master.empty();
  const bool have_slave_meshset = !tie_slave.empty();
 
  if (!have_master_meshset && !have_slave_meshset) {
    MoFEMFunctionReturnHot(0);
  }
 
  auto anyRank = [&](const bool local_value) {
    int global_value = 0;
    const int local_int = static_cast<int>(local_value);
    MPI_Allreduce(&local_int, &global_value, 1, MPI_INT, MPI_MAX,
                  mField.get_comm());
    return global_value != 0;
  };
 
  if (SPACE_DIM != 3) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "TIE prism constraints are supported only in 3D");
  }
 
  if (!tie_data)
    tie_data = boost::make_shared<CommonData>();
 
  tie_data->tieRefLevel = 1;
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-tie_ref_level",
                            &tie_data->tieRefLevel, PETSC_NULLPTR);
  if (tie_data->tieRefLevel < 0) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "-tie_ref_level has to be >= 0");
  }
 
  tie_data->dispFieldName = field_name;
  tie_data->dispFieldBase = base;
  tie_data->dispFieldOrder = order;
 
  if (!have_master_meshset || !have_slave_meshset) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Both TIE_MASTER* and TIE_SLAVE* meshsets are required to "
            "initialise TIE constraints");
  }
  // FIXME: this might not be necessary
  Range local_prisms;
  CHKERR mField.get_moab().get_entities_by_type(0, MBPRISM, local_prisms, true);
  Range local_prism_faces;
  if (!local_prisms.empty()) {
    CHKERR mField.get_moab().get_adjacencies(local_prisms, 2, false,
                                             local_prism_faces,
                                             moab::Interface::UNION);
    local_prism_faces = local_prism_faces.subset_by_type(MBTRI);
  }
 
  auto collectTieFaces = [&](const auto &meshsets,
                             boost::shared_ptr<Range> &face_range) {
    MoFEMFunctionBegin;
    face_range = boost::make_shared<Range>();
    for (const auto &meshset_ptr : meshsets) {
      Range tris;
      // TIE blocksets written to file can expose their surface triangles through
      // nested child sets after reload, so collect them recursively.
      CHKERR mField.get_moab().get_entities_by_type(meshset_ptr->meshset, MBTRI,
                                                    tris, true);
      tris = intersect(tris, local_prism_faces);
      face_range->merge(tris);
    }
    MoFEMFunctionReturn(0);
  };
 
  CHKERR collectTieFaces(tie_master, tie_data->tieMasterFaces);
  CHKERR collectTieFaces(tie_slave, tie_data->tieSlaveFaces);
 
  if (tie_data->tieMasterFaces->empty() != tie_data->tieSlaveFaces->empty()) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Local TIE_MASTER/TIE_SLAVE meshsets are inconsistent");
  }
 
  tie_data->tiePrisms = boost::make_shared<Range>();
  Range master_adjacent_prisms, slave_adjacent_prisms;
  CHKERR mField.get_moab().get_adjacencies(*tie_data->tieSlaveFaces, 3, false,
                                           slave_adjacent_prisms,
                                           moab::Interface::UNION);
  CHKERR mField.get_moab().get_adjacencies(*tie_data->tieMasterFaces, 3, false,
                                           master_adjacent_prisms,
                                           moab::Interface::UNION);
  slave_adjacent_prisms = slave_adjacent_prisms.subset_by_type(MBPRISM);
  master_adjacent_prisms = master_adjacent_prisms.subset_by_type(MBPRISM);
 
  if (!tie_data->tieSlaveFaces->empty() && slave_adjacent_prisms.empty()) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "No prism elements are adjacent to TIE_SLAVE meshsets");
  }
 
  if (slave_adjacent_prisms.size() != master_adjacent_prisms.size()) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Different numbers of prisms are adjacent to TIE_SLAVE and "
            "TIE_MASTER meshsets");
  }
 
  *tie_data->tiePrisms = slave_adjacent_prisms;
 
  const bool have_tie_prisms = anyRank(!tie_data->tiePrisms->empty());
  if (!have_tie_prisms) {
    MOFEM_LOG("WORLD", Sev::warning)
        << "TIE_MASTER/TIE_SLAVE meshsets found but no prism elements are adjacent to them. "
        << "TIE constraints cannot be initialised.";
    MoFEMFunctionReturnHot(0);
  }
 
  if (!tie_data->tiePrisms->empty()) {
    MOFEM_LOG("WORLD", Sev::inform)
        << "Initialising TIE constraints on " << tie_data->tiePrisms->size()
        << " prism elements adjacent to the TIE slave/master surfaces.";
  }
 
  if (!mField.check_field("LAMBDA")) 
    CHKERR mField.add_field("LAMBDA", H1, base, SPACE_DIM);
 
  Range tie_prism_trace_faces = *tie_data->tieSlaveFaces;
  tie_prism_trace_faces.merge(*tie_data->tieMasterFaces);
  Range tie_prism_trace_edges;
  CHKERR mField.get_moab().get_adjacencies(tie_prism_trace_faces, 1, false,
                                           tie_prism_trace_edges,
                                           moab::Interface::UNION);
  Range tie_prism_edges;
  CHKERR mField.get_moab().get_adjacencies(*tie_data->tiePrisms, 1, false,
                                           tie_prism_edges,
                                           moab::Interface::UNION);
  const auto tie_prism_lateral_edges =
      subtract(tie_prism_edges, tie_prism_trace_edges);
 
  CHKERR mField.add_ents_to_field_by_type(*tie_data->tiePrisms, MBPRISM,
                                          field_name);
  CHKERR mField.add_ents_to_field_by_type(*tie_data->tiePrisms, MBPRISM,
                                          "GEOMETRY");
  CHKERR mField.set_field_order(tie_prism_trace_edges, "GEOMETRY", 2);
 
  // LAMBDA is a classical surface multiplier. It lives only on the slave
  // interface, while the prism FE uses that trace space to constrain slave and
  // master displacement traces. Work on the already collected recursive slave
  // face range directly, since reloaded Cubit blocksets can expose triangles
  // only through nested child sets.
  CHKERR mField.add_ents_to_field_by_type(*tie_data->tieSlaveFaces, MBTRI,
                                          "LAMBDA");
  CHKERR comm_interface->synchroniseFieldEntities("LAMBDA", 0);
 
  Range verts;
  CHKERR mField.get_moab().get_entities_by_type(0, MBVERTEX, verts, true);
  Range faces;
  CHKERR mField.get_moab().get_entities_by_type(0, MBTRI, faces, true);
  Range edges;
  CHKERR mField.get_moab().get_adjacencies(faces, 1, true, edges,
                                           moab::Interface::UNION);
  CHKERR mField.add_ents_to_field_by_type(verts, MBVERTEX, "U");
  CHKERR comm_interface->synchroniseEntities(edges);
  CHKERR mField.add_ents_to_field_by_type(edges, MBEDGE, "U");
  CHKERR mField.add_ents_to_field_by_type(faces, MBTRI, "U");
  CHKERR comm_interface->synchroniseFieldEntities("U", 0);
 
 
  CHKERR mField.set_field_order(0, MBVERTEX, field_name, 1);
  CHKERR mField.set_field_order(0, MBEDGE, field_name, order);
  CHKERR mField.set_field_order(0, MBTRI, field_name, order);
  CHKERR mField.set_field_order(tie_prism_lateral_edges, field_name, 1);
  CHKERR mField.set_field_order(0, MBVERTEX, "LAMBDA", 1);
  CHKERR mField.set_field_order(0, MBEDGE, "LAMBDA", order - 1);
  CHKERR mField.set_field_order(0, MBTRI, "LAMBDA", order - 1);
 
  CHKERR mField.build_fields();
 
  if (!mField.check_finite_element(tie_data->tiePrismFeName)) {
    CHKERR mField.add_finite_element(tie_data->tiePrismFeName);
    CHKERR mField.modify_finite_element_add_field_row(tie_data->tiePrismFeName,
                                                      field_name);
    CHKERR mField.modify_finite_element_add_field_col(tie_data->tiePrismFeName,
                                                      field_name);
    CHKERR mField.modify_finite_element_add_field_row(tie_data->tiePrismFeName,
                                                       "LAMBDA");
    CHKERR mField.modify_finite_element_add_field_col(tie_data->tiePrismFeName,
                                                      "LAMBDA");
    CHKERR mField.modify_finite_element_add_field_data(tie_data->tiePrismFeName,
                                                       field_name);
    CHKERR mField.modify_finite_element_add_field_data(tie_data->tiePrismFeName,
                                                       "LAMBDA");
    CHKERR mField.modify_finite_element_add_field_data(tie_data->tiePrismFeName,
                                                       "GEOMETRY");
  }
 
  CHKERR mField.modify_finite_element_add_field_row(simple->getDomainFEName(),
                                                    "LAMBDA");
  CHKERR mField.modify_finite_element_add_field_col(simple->getDomainFEName(),
                                                    "LAMBDA");
  CHKERR mField.modify_finite_element_add_field_data(simple->getDomainFEName(),
                                                     "LAMBDA");
 
  CHKERR mField.add_ents_to_finite_element_by_type(*tie_data->tiePrisms,
                                                   MBPRISM,
                                                   tie_data->tiePrismFeName);
  CHKERR mField.build_finite_elements(simple->getDomainFEName());
  CHKERR mField.build_finite_elements(tie_data->tiePrismFeName);
 
  auto &other_fes = simple->getOtherFiniteElements();
  other_fes.push_back(tie_data->tiePrismFeName);
 
  CHKERR DMMoFEMAddElement(simple->getDM(), tie_data->tiePrismFeName);
  CHKERR simple->reSetUp(true);
 
  tie_data->hasTieConstraints = true;
  MoFEMFunctionReturn(0);
}
 
inline MoFEMErrorCode setupTieConstraintOperators(
    MoFEM::Interface &mField, boost::shared_ptr<CommonData> tie_data) {
  MoFEMFunctionBegin;
 
  if (!tie_data || !tie_data->hasTieConstraints) {
    MoFEMFunctionReturnHot(0);
  }
 
  if (!tie_data->tiePrismRhsFE) {
    tie_data->tiePrismRhsFE =
        boost::make_shared<FlatPrismElementForcesAndSourcesCore>(mField);
    tie_data->tiePrismRhsFE->getRuleHook = [](int, int, int) { return -1; };
    tie_data->tiePrismRhsFE->setRuleHook =
        SetIntegrationOnTiePrismFaces(tie_data->tieRefLevel);
    tie_data->tiePrismRhsFE->getOpPtrVector().push_back(
        new OpCalculateTraceVectorFieldValues<SPACE_DIM>(
            "LAMBDA", tie_data->lambdaAtGaussPts, MBVERTEX));
    tie_data->tiePrismRhsFE->getOpPtrVector().push_back(
        new OpTiePrismRhsU(tie_data->dispFieldName,
                           tie_data->lambdaAtGaussPts));
    tie_data->tiePrismRhsFE->getOpPtrVector().push_back(
        new OpCalculateTraceVectorFieldValues<SPACE_DIM>(
            tie_data->dispFieldName, tie_data->dispAtGaussPts, MBVERTEX));
    tie_data->tiePrismRhsFE->getOpPtrVector().push_back(
        new OpTiePrismRhs(tie_data->dispAtGaussPts));
  }
 
  if (!tie_data->tiePrismLhsFE) {
    tie_data->tiePrismLhsFE =
        boost::make_shared<FlatPrismElementForcesAndSourcesCore>(mField);
    tie_data->tiePrismLhsFE->getRuleHook = [](int, int, int) { return -1; };
    tie_data->tiePrismLhsFE->setRuleHook =
        SetIntegrationOnTiePrismFaces(tie_data->tieRefLevel);
    tie_data->tiePrismLhsFE->getOpPtrVector().push_back(
        new OpTiePrismLhs(tie_data->dispFieldName, "LAMBDA", true));
  }
 
  MoFEMFunctionReturn(0);
}
 
inline MoFEMErrorCode setupTieSolver(MoFEM::Interface &mField,
                                     boost::shared_ptr<CommonData> tie_data,
                                     SmartPetscObj<KSP> solver) {
  MoFEMFunctionBegin;
 
  if (!tie_data || !tie_data->hasTieConstraints)
    MoFEMFunctionReturnHot(0);
 
  CHKERR setupTieConstraintOperators(mField, tie_data);
 
  DM dm;
  CHKERR KSPGetDM(solver, &dm);
  boost::shared_ptr<FEMethod> null;
  CHKERR DMMoFEMKSPSetComputeOperators(dm, tie_data->tiePrismFeName,
                                       tie_data->tiePrismLhsFE, null, null);
  CHKERR DMMoFEMKSPSetComputeRHS(dm, tie_data->tiePrismFeName,
                                 tie_data->tiePrismRhsFE, null, null);
 
  MoFEMFunctionReturn(0);
}
 
inline MoFEMErrorCode setupTieSolver(MoFEM::Interface &mField,
                                     boost::shared_ptr<CommonData> tie_data,
                                     SmartPetscObj<TS> solver) {
  MoFEMFunctionBegin;
 
  if (!tie_data || !tie_data->hasTieConstraints)
    MoFEMFunctionReturnHot(0);
 
  CHKERR setupTieConstraintOperators(mField, tie_data);
 
  DM dm;
  CHKERR TSGetDM(solver, &dm);
  boost::shared_ptr<FEMethod> null;
  CHKERR DMMoFEMTSSetIJacobian(dm, tie_data->tiePrismFeName,
                               tie_data->tiePrismLhsFE, null, null);
  CHKERR DMMoFEMTSSetIFunction(dm, tie_data->tiePrismFeName,
                               tie_data->tiePrismRhsFE, null, null);
 
  MoFEMFunctionReturn(0);
}
 
} // namespace ElasticTie
 
#endif // __ELASTIC_TIE_HPP__