ElasticOperators.cpp

Go to the documentation of this file.

/** \file ElasticOperators.cpp
 */
 
/* This file is part of MoFEM.
 * MoFEM is free software: you can redistribute it and/or modify it under
 * the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 *
 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 * License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
 
#include <MoFEM.hpp>
using namespace MoFEM;
using namespace FTensor;
 
#include <BasicFiniteElements.hpp>
#include <AnalyticalSurfaces.hpp>
#include <RigidBodies.hpp>
#include <BlockMatData.hpp>
#include <MultifieldPlasticity.hpp>
#include <BasicFeTools.hpp>
#include <MatrixFunction.hpp>
 
#include <ElasticOperators.hpp>
#include <chrono> //for debugging
 
struct MeasureTime {
  int time;
  chrono::high_resolution_clock::time_point start;
  chrono::high_resolution_clock::time_point stop;
  MeasureTime() { start = chrono::high_resolution_clock::now(); }
  ~MeasureTime() {
    stop = chrono::high_resolution_clock::now();
    auto duration = chrono::duration_cast<chrono::microseconds>(stop - start);
    cout << "Time taken by function: " << duration.count() << " ms." << endl;
  }
};
 
namespace OpElasticTools {
 
OpStrain::OpStrain(const std::string field_name,
                   boost::shared_ptr<CommonData> common_data_ptr)
    : DomainEleOp(field_name, field_name, DomainEleOp::OPROW),
      commonDataPtr(common_data_ptr) {
  // Operator is only executed for vertices
  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
 
//! [Calculate strain]
MoFEMErrorCode OpStrain::doWork(int side, EntityType type, EntData &data) {
  MoFEMFunctionBegin;
  const size_t nb_gauss_pts = commonDataPtr->mGradPtr->size2();
  commonDataPtr->mStrainPtr->resize(6, nb_gauss_pts);
  auto t_grad = getFTensor2FromMat<3, 3>(*(commonDataPtr->mGradPtr));
  auto t_strain = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStrainPtr));
  commonDataPtr->detF->resize(nb_gauss_pts, false);
  auto t_detF = getFTensor0FromVec(*(commonDataPtr->detF));
 
  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
    t_strain(i, j) = (t_grad(i, j) || t_grad(j, i)) / 2;
    t_detF = 1;
    ++t_detF;
    ++t_grad;
    ++t_strain;
  }
 
  MoFEMFunctionReturn(0);
}
 
OpLogStrain::OpLogStrain(const std::string field_name,
                         boost::shared_ptr<CommonData> common_data_ptr)
    : DomainEleOp(field_name, field_name, DomainEleOp::OPROW),
      commonDataPtr(common_data_ptr) {
  // Operator is only executed for vertices
  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
 
//! [Calculate strain]
MoFEMErrorCode OpLogStrain::doWork(int side, EntityType type, EntData &data) {
  MoFEMFunctionBegin;
  const size_t nb_gauss_pts = commonDataPtr->mGradPtr->size2();
  commonDataPtr->mStrainPtr->resize(6, nb_gauss_pts);
  commonDataPtr->mDeformationGradient->resize(9, nb_gauss_pts, false);
  commonDataPtr->matC->resize(6, nb_gauss_pts, false);
  commonDataPtr->matEigenVal->resize(3, nb_gauss_pts, false);
  commonDataPtr->matEigenVector->resize(9, nb_gauss_pts, false);
  commonDataPtr->detF->resize(nb_gauss_pts, false);
 
  // double dt = 0.;
  // auto get_dt = [&]() {
  //   double dt;
  //   CHKERR TSGetTime(getFEMethod()->ts, &dt);
  //   // CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
  //   dt = dt;
  //   return dt;
  // };
  // if ((getDataCtx() & PetscData::CtxSetTime).any())
  //   dt = get_dt();
 
  auto t_grad = getFTensor2FromMat<3, 3>(*(commonDataPtr->mGradPtr));
  auto F = getFTensor2FromMat<3, 3>(*(commonDataPtr->mDeformationGradient));
  auto C = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->matC));
  auto t_eigen = getFTensor1FromMat<3>(*(commonDataPtr->matEigenVal));
  auto t_eigen_vec = getFTensor2FromMat<3, 3>(*(commonDataPtr->matEigenVector));
  auto t_detF = getFTensor0FromVec(*(commonDataPtr->detF));
  auto t_strain = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStrainPtr));
 
  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
    // auto R = get_rotation_R((*cache).omega, dt);
    // F(i, j) = (t_grad(i, k) + kronecker_delta(i, k)) * R(k, j);
    F(i, j) = t_grad(i, j) + kronecker_delta(i, j);
 
    double det;
    CHKERR determinantTensor3by3(F, det);
    t_detF = det;
 
    // test
    // if (getTStime() > 0.2) {
    //   F(k, i) = 0.;
    //   F(0, 0) = 1;
    //   F(1, 1) = 1;
    //   F(2, 2) = 3;
    //   F(0, 1) = 2;
    //   F(1, 0) = 2;
    // }
 
    C(i, j) = F(k, i) ^ F(k, j);
 
    CHKERR get_eigen_val_and_proj_lapack(C, t_eigen, t_eigen_vec);
    auto lnC = get_ln_X_lapack(C, t_eigen, t_eigen_vec);
    t_strain(i, j) = 0.5 * lnC(i, j);
 
    ++t_eigen;
    ++t_eigen_vec;
    ++C;
    ++t_detF;
    ++F;
    ++t_grad;
    ++t_strain;
  }
 
  MoFEMFunctionReturn(0);
}
OpStress::OpStress(const std::string field_name,
                   boost::shared_ptr<CommonData> common_data_ptr, boost::shared_ptr<MatrixDouble> mat_D_Ptr)
    : DomainEleOp(field_name, field_name, DomainEleOp::OPROW),
      commonDataPtr(common_data_ptr), matDPtr(mat_D_Ptr) {
  // Operator is only executed for vertices
  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
 
//! [Calculate stress]
MoFEMErrorCode OpStress::doWork(int side, EntityType type, EntData &data) {
  MoFEMFunctionBegin;
  const size_t nb_gauss_pts = commonDataPtr->mStrainPtr->size2();
  commonDataPtr->mStressPtr->resize(6, nb_gauss_pts);
  auto t_D = getFTensor4DdgFromMat<3, 3, 0>(*matDPtr);
  auto t_strain = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStrainPtr));
  auto t_stress = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStressPtr));
 
  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
    t_stress(i, j) = t_D(i, j, k, l) * t_strain(k, l);
    ++t_strain;
    ++t_stress;
  }
 
  MoFEMFunctionReturn(0);
}
template <bool TANGENT>
OpCalculateLogStressTangent<TANGENT>::OpCalculateLogStressTangent(
    const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr,
    boost::shared_ptr<MatrixDouble> mat_D_ptr)
    : DomainEleOp(field_name, field_name, DomainEleOp::OPROW),
      commonDataPtr(common_data_ptr), matDPtr(mat_D_ptr) {
  // Operator is only executed for vertices
  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
 
template <bool TANGENT>
MoFEMErrorCode OpCalculateLogStressTangent<TANGENT>::doWork(int side,
                                                            EntityType type,
                                                            EntData &data) {
  MoFEMFunctionBegin;
 
  const size_t nb_gauss_pts = data.getN().size1();
  auto t_D = getFTensor4DdgFromMat<3, 3, 0>(*matDPtr);
  auto t_stress = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStressPtr));
  commonDataPtr->mPiolaStressPtr->resize(9, nb_gauss_pts, false);
  auto t_piola = getFTensor2FromMat<3, 3>(*(commonDataPtr->mPiolaStressPtr));
  auto F = getFTensor2FromMat<3, 3>(*(commonDataPtr->mDeformationGradient));
  auto C = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->matC));
  auto t_eigen_vec = getFTensor2FromMat<3, 3>(*(commonDataPtr->matEigenVector));
  auto t_eigen = getFTensor1FromMat<3>(*(commonDataPtr->matEigenVal));
 
  MatrixDouble &dP_dF = *(commonDataPtr->materialTangent);
  dP_dF.resize(81, nb_gauss_pts, false);
  auto t_dP_dF = getFTensor4FromMat<3, 3, 3, 3>(dP_dF);
 
  MatrixDouble &dE_dF = *(commonDataPtr->dE_dF_mat);
  dE_dF.resize(81, nb_gauss_pts, false);
  auto t_dE_dF = getFTensor4FromMat<3, 3, 3, 3>(dE_dF);
 
  constexpr auto t_kd = Kronecker_Delta<int>();
  Ddg<double, 3, 3> TL;
  // Ddg<double, 3, 3> TL_test;
  Tensor2_symmetric<double, 3> S;
 
  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
    FTensor::Tensor4<double, 3, 3, 3, 3> dC_dF;
    dC_dF(i, j, k, l) = (t_kd(i, l) * F(k, j)) + (t_kd(j, l) * F(k, i));
    auto dlogC_dC = calculate_lagrangian_moduli_TL<TANGENT>(t_stress, t_eigen,
                                                            t_eigen_vec, TL);
 
    t_dE_dF(i, j, k, l) = dlogC_dC(i, j, m, n) * dC_dF(m, n, k, l);
    t_dE_dF(i, j, k, l) *= 0.5;
 
    if (TANGENT) {
      S(k, l) = t_stress(i, j) * dlogC_dC(i, j, k, l);
      t_piola(i, l) = F(i, k) * S(k, l);
 
      FTensor::Ddg<double, 3, 3> P_D_P_plus_TL;
      P_D_P_plus_TL(i, j, k, l) =
          TL(i, j, k, l) +
          (dlogC_dC(i, j, o, p) * t_D(o, p, m, n)) * dlogC_dC(m, n, k, l);
      P_D_P_plus_TL(i, j, k, l) *= 0.5;
 
      t_dP_dF(i, j, m, n) = t_kd(i, m) * (t_kd(k, n) * S(k, j));
      t_dP_dF(i, j, m, n) +=
          F(i, k) * (P_D_P_plus_TL(k, j, o, p) * dC_dF(o, p, m, n));
 
      // auto Piola = to_non_symm(t_stress);
      // auto D1_from_fddsa321p = getDTensorFunction(*commonDataPtr, F, Piola);
      // auto test_lukasz =
      //     test_projection_lukasz(C, t_stress, t_eigen, t_eigen_vec, TL_test);
      // test_lukasz(i, j, k, l) *= 2;
      // TL_test(i, j, k, l) *= 4;
      // string wait;
    } else {
      // auto dlogC_dC = calculate_lagrangian_moduli_TL<false>(
      //     t_stress, t_eigen, t_eigen_vec, TL);
      S(k, l) = t_stress(i, j) * dlogC_dC(i, j, k, l);
      t_piola(i, l) = F(i, k) * S(k, l);
    }
    if (false && TANGENT) {
      if (getTStime() > 0.2) {
        {
          Ddg<double, 3, 3> TL3_lukasz;
          cout << "NEW IMPLEMENTATION" << endl;
          MeasureTime time1;
          for (int ii = 0; ii != 1000; ii++)
            auto P_proj_lukasz = test_projection_lukasz(
                C, t_stress, t_eigen, t_eigen_vec, TL3_lukasz);
        }
        Ddg<double, 3, 3> dummy33;
        {
          cout << "NEW IMPLEMENTATION" << endl;
          MeasureTime time1;
          for (int ii = 0; ii != 1000; ii++)
            auto P_proj_ = calculate_lagrangian_moduli_TL<true>(
                t_stress, t_eigen, t_eigen_vec, dummy33);
        }
        {
          cout << "(OLD) FINITE DIFFERENCE" << endl;
          auto Piola = to_non_symm(t_stress);
          MeasureTime time1;
          for (int ii = 0; ii != 1000; ii++)
            auto D1_from_fddsap = getDTensorFunction(*commonDataPtr, F, Piola);
        }
        {
          cout << "EIGENVALUES (LAPACK)" << endl;
          MeasureTime time1;
          for (int ii = 0; ii != 1000; ii++)
            CHKERR get_eigen_val_and_proj_lapack(C, t_eigen, t_eigen_vec);
        }
      }
 
      // double sum = 0;
      // double sum2 = 0;
      // double sum3 = 0;
      // double sum0 = 0;
      // for (int i = 0; i != 3; ++i)
      //   for (int j = 0; j != 3; ++j)
      //     for (int k = 0; k != 3; ++k)
      //       for (int l = 0; l != 3; ++l) {
      //         // sum0 += abs(P3_prop_lag(i, j, k, l) - dlogC_dC_test(i, j, k,
      //         l));
      //         // sum += abs(P_Proj_lukasz(i, j, k, l) - dlogC_dC_test(i, j,
      //         k, l)); sum2 += abs(dP_dF(i, j, k, l) - D1(i, j, k, l)); sum3
      //         +=
      //             abs(myDFin_lukasz(i, j, k, l) - D1_from_fd_LAG(i, j, k,
      //             l));
      //         // if (abs(myDFinFinFin(i, j, k, l) - D1(i, j, k, l)) > 1e-6)
      //         //   cout << i << j << k << l << " "
      //         //        << myDFinFinFin(i, j, k, l) / D1(i, j, k, l) << endl;
      //       }
      if (getTStime() > 0.2)
        string wait;
    }
 
    ++t_piola;
    ++t_eigen;
    ++t_eigen_vec;
    ++t_dP_dF;
    ++t_dE_dF;
    ++t_stress;
    ++F;
    ++C;
  }
 
  MoFEMFunctionReturn(0);
}
 
OpCalculateNeoHookeStressTangent::OpCalculateNeoHookeStressTangent(
    const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
    : DomainEleOp(field_name, field_name, DomainEleOp::OPROW),
      commonDataPtr(common_data_ptr) {
  // Operator is only executed for vertices
  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
MoFEMErrorCode OpCalculateNeoHookeStressTangent::doWork(int side,
                                                        EntityType type,
                                                        EntData &data) {
  MoFEMFunctionBegin;
 
  const size_t nb_gauss_pts = data.getN().size1();
  auto t_stress = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStressPtr));
  commonDataPtr->mPiolaStressPtr->resize(9, nb_gauss_pts, false);
  auto t_piola = getFTensor2FromMat<3, 3>(*(commonDataPtr->mPiolaStressPtr));
  auto F = getFTensor2FromMat<3, 3>(*(commonDataPtr->mDeformationGradient));
  auto t_detF = getFTensor0FromVec(*(commonDataPtr->detF));
 
  MatrixDouble &dP_dF = *(commonDataPtr->materialTangent);
  dP_dF.resize(81, nb_gauss_pts, false);
  auto t_dP_dF = getFTensor4FromMat<3, 3, 3, 3>(dP_dF);
 
  MatrixDouble &dE_dF = *(commonDataPtr->dE_dF_mat);
  dE_dF.resize(81, nb_gauss_pts, false);
  auto t_dE_dF = getFTensor4FromMat<3, 3, 3, 3>(dE_dF);
  const bool is_lhs = (int)getTSCtx() == 3;
  constexpr auto t_kd = Kronecker_Delta<int>();
 
  const double lambda = LAMBDA((*cache).young_modulus, (*cache).poisson);
  const double mu = MU((*cache).young_modulus, (*cache).poisson);
 
  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
 
    Tensor2<double, 3, 3> invF;
    CHKERR invertTensor3by3(F, t_detF, invF);
    const double log_det = log(t_detF);
 
    if (is_lhs) {
      // compute tangent
      const double coef = mu - lambda * log_det;
      t_dP_dF(i, J, k, L) =
          invF(J, i) * (invF(L, k) * lambda) + invF(L, i) * (invF(J, k) * coef);
 
      t_dP_dF(0, 0, 0, 0) += mu;
      t_dP_dF(0, 1, 0, 1) += mu;
      t_dP_dF(0, 2, 0, 2) += mu;
      t_dP_dF(1, 0, 1, 0) += mu;
      t_dP_dF(1, 1, 1, 1) += mu;
      t_dP_dF(1, 2, 1, 2) += mu;
      t_dP_dF(2, 0, 2, 0) += mu;
      t_dP_dF(2, 1, 2, 1) += mu;
      t_dP_dF(2, 2, 2, 2) += mu;
 
    } else {
 
      // compute only stress
      const double var = lambda * log_det - mu;
      t_piola(i, J) = mu * F(i, J) + var * invF(J, i);
    }
 
    ++t_detF;
    ++t_piola;
    ++t_dP_dF;
    ++t_stress;
    ++F;
  }
 
  MoFEMFunctionReturn(0);
}
 
//! [Calculate stress]
template <bool LOGSTRAIN, bool REACTIONS>
OpInternalForceRhs<LOGSTRAIN, REACTIONS>::OpInternalForceRhs(
    const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
    : DomainEleOpAssembly(field_name, field_name, DomainEleOp::OPROW),
      commonDataPtr(common_data_ptr) {}
 
//! [Internal force]
template <bool LOGSTRAIN, bool REACTIONS>
MoFEMErrorCode
OpInternalForceRhs<LOGSTRAIN, REACTIONS>::iNtegrate(EntData &data) {
  MoFEMFunctionBegin;
 
  const size_t nb_dofs = data.getIndices().size();
  if (nb_dofs) {
 
    const size_t nb_base_functions = data.getN().size2();
    if (3 * nb_base_functions < nb_dofs)
      SETERRQ(
          PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
          "Number of DOFs is larger than number of base functions on entity");
 
    const size_t nb_gauss_pts = data.getN().size1();
 
    auto t_w = getFTensor0IntegrationWeight();
    auto t_stress =
        getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStressPtr));
    auto t_piola = getFTensor2FromMat<3, 3>(*(commonDataPtr->mPiolaStressPtr));
    auto t_diff_base = data.getFTensor1DiffN<3>();
 
    for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
 
      double alpha = getMeasure() * t_w;
      auto t_nf = getFTensor1FromArray<3, 3>(locF);
      auto out_stress = to_non_symm(t_stress);
 
      if (LOGSTRAIN) {
        out_stress(i, j) = t_piola(i, j);
        ++t_piola;
      }
 
      size_t bb = 0;
      for (; bb != nb_dofs / 3; ++bb) {
        t_nf(i) += alpha * t_diff_base(j) * out_stress(i, j);
        ++t_diff_base;
        ++t_nf;
      }
      for (; bb < nb_base_functions; ++bb)
        ++t_diff_base;
 
      ++t_stress;
      ++t_w;
    }
  }
  MoFEMFunctionReturn(0);
}
 
template <bool LOGSTRAIN, bool REACTIONS>
MoFEMErrorCode
OpInternalForceRhs<LOGSTRAIN, REACTIONS>::aSsemble(EntData &data) {
  MoFEMFunctionBeginHot;
  if (REACTIONS) {
    auto react_ents = commonDataPtr->reactionEnts;
    EntityHandle ent = getFEEntityHandle();
 
    auto &vec = commonDataPtr->reactionsVec;
    double *react_ptr;
    CHKERR VecGetArray(vec, &react_ptr);
 
    int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
    int ii = 0;
    for (auto &dof : data.getFieldDofs()) {
      auto ent = dof->getEnt();
      double tag_val[rank];
      auto idx = dof->getDofCoeffIdx();
      auto &m_field = getPtrFE()->mField;
      CHKERR m_field.get_moab().tag_get_data(commonDataPtr->reactionTag, &ent,
                                             1, tag_val);
      if (dof->getEntType() == MBVERTEX) {
        const double react = this->locF(ii) / (*cache).young_modulus_inv;
        tag_val[idx] += react;
        if (react_ents.find(ent) != react_ents.end())
          react_ptr[idx] += react;
      }
      ii++;
      CHKERR m_field.get_moab().tag_set_data(commonDataPtr->reactionTag, &ent,
                                             1, &tag_val);
    }
    CHKERR VecRestoreArray(vec, &react_ptr);
  } else
    CHKERR VecSetValues<EssentialBcStorage>(
        this->getKSPf(), data, &*this->locF.data().begin(), ADD_VALUES);
 
  MoFEMFunctionReturnHot(0);
}
 
OpSetMaterialBlock::OpSetMaterialBlock(
    const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr,
    std::map<EntityHandle, int> &map_block_tets)
    : DomainEleOp(field_name, field_name, DomainEleOp::OPROW),
      commonDataPtr(common_data_ptr), mapBlockTets(map_block_tets) {
  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
 
MoFEMErrorCode OpSetMaterialBlock::doWork(int side, EntityType type,
                                          EntData &data) {
  MoFEMFunctionBegin;
  if (mapBlockTets.empty())
    MoFEMFunctionReturnHot(0);
 
  auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
  const auto id = mapBlockTets.at(fe_ent);
  cache = &mat_blocks.at(id);
  *commonDataPtr->mtD = (*cache).mtD;
 
  // TODO:
  // FIXME: include params for plasticity
  MoFEMFunctionReturn(0);
}
 
OpSetMaterialBlockBoundary::OpSetMaterialBlockBoundary(
    const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr,
    std::map<EntityHandle, int> &map_block_tets)
    : BoundaryEleOp(field_name, field_name, BoundaryEleOp::OPROW),
      commonDataPtr(common_data_ptr), mapBlockTets(map_block_tets) {
  std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
  doEntities[MBTRI] = doEntities[MBQUAD] = true;
}
 
MoFEMErrorCode OpSetMaterialBlockBoundary::doWork(int side, EntityType type,
                                                  EntData &data) {
  MoFEMFunctionBegin;
  if (mapBlockTets.empty())
    MoFEMFunctionReturnHot(0);
  std::vector<EntityHandle> adj_ents;
  auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
  auto &moab = getPtrFE()->mField.get_moab();
  CHKERR moab.get_adjacencies(&fe_ent, 1, 3, false, adj_ents);
 
  const auto id = mapBlockTets.at(adj_ents[0]);
  cache = &mat_blocks.at(id);
  *commonDataPtr->mtD = (*cache).mtD;
 
  MoFEMFunctionReturn(0);
}
 
OpSaveReactionForces::OpSaveReactionForces(
    MoFEM::Interface &m_field, const std::string field_name,
    moab::Interface &post_proc_mesh, std::vector<EntityHandle> &map_gauss_pts,
    boost::shared_ptr<CommonData> common_data_ptr)
    : DomainEleOp(field_name, field_name, DomainEleOp::OPROW), mField(m_field),
      postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
      commonDataPtr(common_data_ptr) {
  // Operator is only executed for vertices
  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
 
MoFEMErrorCode OpSaveReactionForces::doWork(int side, EntityType type,
                                            EntData &data) {
  MoFEMFunctionBegin;
 
  auto get_tag = [&](const std::string name, size_t size) {
    std::array<double, 9> def;
    std::fill(def.begin(), def.end(), 0);
    Tag th;
    CHKERR postProcMesh.tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE, th,
                                       MB_TAG_CREAT | MB_TAG_SPARSE,
                                       def.data());
    return th;
  };
 
  auto th_react = get_tag("REACTIONS", 3);
  size_t nb_gauss_pts = data.getN().size1();
  VectorDouble vAlues(data.getFieldData().size());
 
  const size_t rank = data.getFieldDofs()[0]->getNbOfCoeffs();
 
  int dd = 0;
  for (auto &dof : data.getFieldDofs()) {
    auto ent = dof->getEnt();
    // ents_vec[dd] = ent;
    const size_t idx = dof->getDofCoeffIdx();
 
    double tag_val[rank];
    CHKERR mField.get_moab().tag_get_data(commonDataPtr->reactionTag, &ent, 1,
                                          tag_val);
    vAlues(dd++) = tag_val[idx];
  }
 
  const void *tags_ptr[mapGaussPts.size()];
  CHKERR postProcMesh.tag_get_by_ptr(th_react, &mapGaussPts[0],
                                     mapGaussPts.size(), tags_ptr);
  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
    for (int rr = 0; rr != rank; rr++) {
      // its ugly, but tested in PostProcOnRefMesh
      const double val = cblas_ddot((vAlues.size() / rank), &(data.getN(gg)[0]),
                                    1, &vAlues(rr), rank);
      ((double *)tags_ptr[gg])[rr] += val;
    }
  }
 
  MoFEMFunctionReturn(0);
}
template <bool LOGSTRAIN>
OpPostProcElastic<LOGSTRAIN>::OpPostProcElastic(
    const std::string field_name, moab::Interface &post_proc_mesh,
    std::vector<EntityHandle> &map_gauss_pts,
    boost::shared_ptr<CommonData> common_data_ptr)
    : DomainEleOp(field_name, field_name, DomainEleOp::OPROW),
      postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
      commonDataPtr(common_data_ptr) {
  // Operator is only executed for vertices
  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
 
//! [Postprocessing]
template <bool LOGSTRAIN>
MoFEMErrorCode OpPostProcElastic<LOGSTRAIN>::doWork(int side, EntityType type,
                                                    EntData &data) {
  MoFEMFunctionBegin;
 
  auto get_tag = [&](const std::string name, size_t size = 9) {
    std::array<double, 9> def;
    std::fill(def.begin(), def.end(), 0);
    Tag th;
    CHKERR postProcMesh.tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE, th,
                                       MB_TAG_CREAT | MB_TAG_SPARSE,
                                       def.data());
    return th;
  };
 
  MatrixDouble3by3 mat(3, 3);
 
  auto set_matrix = [&](auto &t) -> MatrixDouble3by3 & {
    mat.clear();
    for (size_t r = 0; r != 3; ++r)
      for (size_t c = 0; c != 3; ++c)
        mat(r, c) = t(r, c);
    return mat;
  };
 
  auto set_scalar = [&](auto t) -> MatrixDouble3by3 & {
    mat.clear();
    mat(0, 0) = t;
    return mat;
  };
 
  auto set_tag = [&](auto th, auto gg, MatrixDouble3by3 &mat) {
    return postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
                                     &*mat.data().begin());
  };
 
  string strain_tag_name = "STRAIN";
  string stress_tag_name = "CAUCHY_STRESS";
  if (LOGSTRAIN) {
    strain_tag_name = "LOG_STRAIN";
    // stress_tag_name = "PIOLA_KIRCHHOFF_STRESS";
  }
 
  auto th_strain = get_tag(strain_tag_name);
  Tensor2<double, 3, 3> F;
  auto th_grad = get_tag("GRAD");
  auto th_stress = get_tag(stress_tag_name);
  auto th_mises = get_tag("MISES", 1);
 
  size_t nb_gauss_pts = data.getN().size1();
  auto t_grad = getFTensor2FromMat<3, 3>(*(commonDataPtr->mGradPtr));
  auto t_strain = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStrainPtr));
  auto t_stress = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStressPtr));
  auto t_piola = getFTensor2FromMat<3, 3>(*(commonDataPtr->mPiolaStressPtr));
 
  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
 
    t_stress(i, j) /= (*cache).young_modulus_inv;
 
    if (LOGSTRAIN) {
      // auto o_stress = to_non_symm(t_stress);
      t_piola(i, j) /= (*cache).young_modulus_inv;
      F(i, j) = t_grad(i, j) + kronecker_delta(i, j);
      double detF;
      CHKERR determinantTensor3by3(F, detF);
      // cauchy
      t_stress(i, j) = (1. / detF) * (t_piola(i, k) ^ F(j, k));
      ++t_piola;
    }
 
    const double mises =
        sqrt(0.5 * (pow(t_stress(0, 0) - t_stress(1, 1), 2) +
                    pow(t_stress(1, 1) - t_stress(2, 2), 2) +
                    pow(t_stress(2, 2) - t_stress(0, 0), 2) +
                    6 * (pow(t_stress(0, 1), 2) + pow(t_stress(1, 2), 2) +
                         pow(t_stress(2, 0), 2))));
 
    CHKERR set_tag(th_mises, gg, set_scalar(mises));
    // FIXME: add option for more postprocessing
    // CHKERR set_tag(th_grad, gg, set_matrix(t_grad));
    if (!commonDataPtr->isNeoHook)
      CHKERR set_tag(th_strain, gg, set_matrix(t_strain));
    CHKERR set_tag(th_stress, gg, set_matrix(t_stress));
 
    ++t_grad;
    ++t_strain;
    ++t_stress;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSchurBeginBoundary::doWork(int row_side, EntityType row_type,
                                    EntData &data) {
  MoFEMFunctionBegin;
  if (row_type != MBTRI && row_type != MBQUAD)
    MoFEMFunctionReturnHot(0);
  
  // MoFEMFunctionReturnHot(0); //FIXME:
  if (!commonDataPtr->blockMatContainerPtr)
    MoFEMFunctionReturnHot(0);
  auto &bmc = *commonDataPtr->blockMatContainerPtr;
  for (auto &bit : bmc)
    bit.unSetAtElement();
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSchurEndBoundary::doWork(int row_side, EntityType row_type,
                                  EntData &data) {
  MoFEMFunctionBegin;
  if (row_type != MBTRI && row_type != MBQUAD)
    MoFEMFunctionReturnHot(0);
  
  // MoFEMFunctionReturnHot(0); //FIXME:
  
  if (commonDataPtr->blockMatContainerPtr) {
 
    auto SEpEp = commonDataPtr->SEpEp;
    auto aoSEpEp = commonDataPtr->aoSEpEp;
    auto STauTau = commonDataPtr->STauTau;
    auto aoSTauTau = commonDataPtr->aoSTauTau;
 
    if (SEpEp) {
 
      auto find_block = [&](BlockMatContainer &add_bmc, auto &row_name,
                            auto &col_name, auto row_side, auto col_side,
                            auto row_type, auto col_type) {
        return add_bmc.get<0>().find(boost::make_tuple(
            row_name, col_name, row_type, col_type, row_side, col_side));
      };
 
      auto set_block = [&](BlockMatContainer &add_bmc, auto &row_name,
                           auto &col_name, auto row_side, auto col_side,
                           auto row_type, auto col_type, const auto &m,
                           const auto &row_ind, const auto &col_ind) {
        auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
                             row_type, col_type);
        if (it != add_bmc.get<0>().end()) {
          it->setMat(m);
          it->setInd(row_ind, col_ind);
          it->setSetAtElement();
          return it;
        } else {
          auto p_it = add_bmc.insert(BlockMatData(row_name, col_name, row_type,
                                                  col_type, row_side, col_side,
                                                  m, row_ind, col_ind));
          return p_it.first;
        }
      };
 
      auto add_block = [&](BlockMatContainer &add_bmc, auto &row_name,
                           auto &col_name, auto row_side, auto col_side,
                           auto row_type, auto col_type, const auto &m,
                           const auto &row_ind, const auto &col_ind) {
        auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
                             row_type, col_type);
        if (it != add_bmc.get<0>().end()) {
          it->addMat(m);
          it->setInd(row_ind, col_ind);
          it->setSetAtElement();
          return it;
        } else {
          auto p_it = add_bmc.insert(BlockMatData(row_name, col_name, row_type,
                                                  col_type, row_side, col_side,
                                                  m, row_ind, col_ind));
          return p_it.first;
        }
      };
 
      auto assemble_block = [&](auto &bit, Mat S) {
        MoFEMFunctionBeginHot;
        const VectorInt &rind = bit.rowInd;
        const VectorInt &cind = bit.colInd;
        const MatrixDouble &m = bit.M;
 
        CHKERR MatSetValues(S, rind.size(), &*rind.begin(), cind.size(),
                            &*cind.begin(), &*m.data().begin(), ADD_VALUES);
 
        MoFEMFunctionReturnHot(0);
      };
 
      auto create_block_schur_bc = [&](BlockMatContainer &bmc,
                                    BlockMatContainer &add_bmc,
                                    std::string field, AO ao) {
        MoFEMFunctionBeginHot;
        for (auto &bit : add_bmc) {
          bit.unSetAtElement();
          bit.clearMat();
        }
 
        for (auto &bit : bmc) {
          if (bit.setAtElement && bit.rowField != field &&
              bit.colField != field) {
            VectorInt rind = bit.rowInd;
            VectorInt cind = bit.colInd;
            const MatrixDouble &m = bit.M;
 
            // for (int i = 0; i != m.size2() * m.size1(); i++)
            //   if (abs(m.data()[i]) > 1e-8) {
            //     string wait;
            //     cout << "what" << endl;
            //   }
 
            if (ao) {
              CHKERR AOApplicationToPetsc(ao, rind.size(), &*rind.begin());
              CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
            }
            auto it =
                set_block(add_bmc, bit.rowField, bit.colField, bit.rowSide,
                          bit.colSide, bit.rowType, bit.colType, m, rind, cind);
          }
        }
 
        // for (auto &bit_col : bmc) {
        //   if (bit_col.setAtElement && bit_col.rowField == field &&
        //       bit_col.colField != field) {
        //     const MatrixDouble &cm = bit_col.M;
        //     VectorInt cind = bit_col.colInd;
        //     // invMat.resize(inv.size1(), cm.size2(), false);
        //     // noalias(invMat) = prod(inv, cm);
        //     noalias(invMat) = cm;
        //     if (ao)
        //       CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
        //     for (auto &bit_row : bmc) {
        //       if (bit_row.setAtElement && bit_row.rowField != field &&
        //           bit_row.colField == field) {
        //         const MatrixDouble &rm = bit_row.M;
        //         VectorInt rind = bit_row.rowInd;
        //         K.resize(rind.size(), cind.size(), false);
        //         // noalias(K) = prod(rm, invMat);
        //         noalias(K) = rm;
        //         if (ao)
        //           CHKERR AOApplicationToPetsc(ao, rind.size(), &*rind.begin());
        //         auto it =
        //             add_block(add_bmc, bit_row.rowField, bit_col.colField,
        //                       bit_row.rowSide, bit_col.colSide, bit_row.rowType,
        //                       bit_col.colType, K, rind, cind);
        //       }
        //     }
        //   }
        // }
 
        MoFEMFunctionReturnHot(0);
      };
 
      auto assemble_schur_bc = [&](BlockMatContainer &add_bmc, Mat S,
                                bool debug = false) {
        MoFEMFunctionBeginHot;
        for (auto &bit : add_bmc) {
          if (bit.setAtElement)
            CHKERR assemble_block(bit, S);
        }
        if (debug) {
          for (auto &bit : add_bmc) {
            if (bit.setAtElement) {
              std::cerr << "assemble: " << bit.rowField << " " << bit.colField
                        << endl;
              // std::cerr << bit.M << endl;
            }
          }
          std::cerr << std::endl;
        }
        MoFEMFunctionReturnHot(0);
      };
 
       auto precondition_schur = [&](BlockMatContainer &bmc,
                                    BlockMatContainer &add_bmc,
                                    const std::string field,
                                    const MatrixDouble &diag_mat,
                                    const double eps) {
        MoFEMFunctionBeginHot;
 
        for (auto &bit : add_bmc) {
          bit.unSetAtElement();
          bit.clearMat();
        }
 
        for (auto &bit : bmc) {
          if (bit.setAtElement) {
            if (bit.rowField != field || bit.colField != field)
              auto it = set_block(add_bmc, bit.rowField, bit.colField,
                                  bit.rowSide, bit.colSide, bit.rowType,
                                  bit.colType, bit.M, bit.rowInd, bit.colInd);
          }
        }
 
        auto bit =
            bmc.get<1>().find(boost::make_tuple(field, field, row_type, row_type));
        if (bit->setAtElement && bit != bmc.get<1>().end()) {
          auto it = set_block(add_bmc, bit->rowField, bit->colField,
                              bit->rowSide, bit->colSide, bit->rowType,
                              bit->colType, bit->M, bit->rowInd, bit->colInd);
          MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
          m += eps * diag_mat;
        } else {
          auto row_it = bmc.get<3>().lower_bound(field);
          for (; row_it != bmc.get<3>().upper_bound(field); ++row_it) {
            if (row_it->setAtElement) {
              auto it = set_block(add_bmc, field, field, 0, 0, row_type, row_type,
                                  diag_mat, row_it->rowInd, row_it->rowInd);
              MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
              m *= eps;
              break;
            }
          }
          if (row_it == bmc.get<3>().end())
            SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
                     "row field not found %s", field.c_str());
        }
 
        MoFEMFunctionReturnHot(0);
      };
 
      const bool debug = false;
      auto &mat_cont = *commonDataPtr->blockMatContainerPtr;
      auto &mass_mat_tensor = commonDataPtr->massMatTensor;
 
      if (SEpEp) {
        CHKERR create_block_schur_bc(mat_cont, blockMat["BC"], "EP", aoSEpEp);
        // if (getPtrFE()->mField.check_field("SIGMA") && false) {
        //   constexpr double eps = 1e-8;
        //   CHKERR precondition_schur(blockMat["BC"], blockMat["precBC"],
        //   "SIGMA",
        //                             mass_mat_tensor, eps);
        //   CHKERR assemble_schur_bc(blockMat["precBC"], SEpEp, debug);
        //   blockMat["BC"] = blockMat["precBC"];
        //   // CHKERR assemble_schur_bc(blockMat["BC"], SEpEp, debug);
        // } else
        CHKERR assemble_schur_bc(blockMat["BC"], SEpEp, debug);
 
        if (STauTau) {
          CHKERR create_block_schur_bc(blockMat["BC"], blockMat["TAUTAU"],
                                       "TAU", aoSTauTau);
          CHKERR assemble_schur_bc(blockMat["TAUTAU"], STauTau, debug);
        }
      }
 
      blockMat.clear();
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSchurBegin::doWork(int row_side, EntityType row_type,
                                    EntData &data) {
  MoFEMFunctionBegin;
  if (row_type != MBTET && row_type != MBHEX)
    MoFEMFunctionReturnHot(0);
  if (!commonDataPtr->blockMatContainerPtr)
    MoFEMFunctionReturnHot(0);
  auto &bmc = *commonDataPtr->blockMatContainerPtr;
  for (auto &bit : bmc)
    bit.unSetAtElement();
  // bmc.clear();
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSchurEnd::doWork(int row_side, EntityType row_type,
                                  EntData &data) {
  MoFEMFunctionBegin;
  if (row_type != MBTET && row_type != MBHEX)
    MoFEMFunctionReturnHot(0);
  if (commonDataPtr->blockMatContainerPtr) {
 
    auto SEpEp = commonDataPtr->SEpEp;
    auto aoSEpEp = commonDataPtr->aoSEpEp;
    auto STauTau = commonDataPtr->STauTau;
    auto aoSTauTau = commonDataPtr->aoSTauTau;
 
    if (SEpEp) {
 
      auto find_block = [&](BlockMatContainer &add_bmc, auto &row_name,
                            auto &col_name, auto row_side, auto col_side,
                            auto row_type, auto col_type) {
        return add_bmc.get<0>().find(boost::make_tuple(
            row_name, col_name, row_type, col_type, row_side, col_side));
      };
 
      auto set_block = [&](BlockMatContainer &add_bmc, auto &row_name,
                           auto &col_name, auto row_side, auto col_side,
                           auto row_type, auto col_type, const auto &m,
                           const auto &row_ind, const auto &col_ind) {
        auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
                             row_type, col_type);
        if (it != add_bmc.get<0>().end()) {
          it->setMat(m);
          it->setInd(row_ind, col_ind);
          it->setSetAtElement();
          return it;
        } else {
          auto p_it = add_bmc.insert(BlockMatData(row_name, col_name, row_type,
                                                  col_type, row_side, col_side,
                                                  m, row_ind, col_ind));
          return p_it.first;
        }
      };
 
      auto add_block = [&](BlockMatContainer &add_bmc, auto &row_name,
                           auto &col_name, auto row_side, auto col_side,
                           auto row_type, auto col_type, const auto &m,
                           const auto &row_ind, const auto &col_ind) {
        auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
                             row_type, col_type);
        if (it != add_bmc.get<0>().end()) {
          it->addMat(m);
          it->setInd(row_ind, col_ind);
          it->setSetAtElement();
          return it;
        } else {
          auto p_it = add_bmc.insert(BlockMatData(row_name, col_name, row_type,
                                                  col_type, row_side, col_side,
                                                  m, row_ind, col_ind));
          return p_it.first;
        }
      };
 
      auto assemble_block = [&](auto &bit, Mat S) {
        MoFEMFunctionBeginHot;
        const VectorInt &rind = bit.rowInd;
        const VectorInt &cind = bit.colInd;
        const MatrixDouble &m = bit.M;
 
        CHKERR MatSetValues(S, rind.size(), &*rind.begin(), cind.size(),
                            &*cind.begin(), &*m.data().begin(), ADD_VALUES);
 
        MoFEMFunctionReturnHot(0);
      };
 
      auto invert_symm_mat = [&](MatrixDouble &m, auto &inv) {
        MoFEMFunctionBeginHot;
        const int nb = m.size1();
 
        inv.resize(nb, nb, false);
        inv.clear();
        for (int cc = 0; cc != nb; ++cc)
          inv(cc, cc) = -1;
 
        iPIV.resize(nb, false);
        lapackWork.resize(nb * nb, false);
        const auto info = lapack_dsysv('L', nb, nb, &*m.data().begin(), nb,
                                       &*iPIV.begin(), &*inv.data().begin(), nb,
                                       &*lapackWork.begin(), nb * nb);
        // if (info != 0)
        //   SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
        //            "Can not invert matrix info = %d", info);
 
        MoFEMFunctionReturnHot(0);
      };
 
      auto invert_symm_schur = [&](BlockMatContainer &bmc, std::string field,
                                   auto &inv) {
        MoFEMFunctionBeginHot;
 
        auto bit =
            bmc.get<1>().find(boost::make_tuple(field, field, row_type, row_type));
        if (bit != bmc.get<1>().end()) {
 
          auto &m = *const_cast<MatrixDouble *>(&(bit->M));
          CHKERR invert_symm_mat(m, inv);
 
        } else
          SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
                   "%s matrix not found", field.c_str());
 
        MoFEMFunctionReturnHot(0);
      };
 
      // auto invert_nonsymm_mat = [&](MatrixDouble &m, auto &inv) {
      //   MoFEMFunctionBeginHot;
 
      //   const int nb = m.size1();
 
      //   MatrixDouble trans_m = trans(m);
      //   MatrixDouble trans_inv;
      //   trans_inv.resize(nb, nb, false);
      //   trans_inv.clear();
      //   for (int c = 0; c != nb; ++c)
      //     trans_inv(c, c) = -1;
 
      //   iPIV.resize(nb, false);
      //   const auto info =
      //       lapack_dgesv(nb, nb, &*trans_m.data().begin(), nb,
      //       &*iPIV.begin(),
      //                    &*trans_inv.data().begin(), nb);
      //   if (info != 0)
      //     SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
      //              "Can not invert matrix info = %d", info);
 
      //   inv.resize(nb, nb, false);
      //   noalias(inv) = trans(trans_inv);
 
      //   MoFEMFunctionReturnHot(0);
      // };
 
      // auto invert_nonsymm_schur = [&](BlockMatContainer &bmc,
      //                                 std::string field, auto &inv,
      //                                 const bool debug = false) {
      //   MoFEMFunctionBeginHot;
 
      //   auto bit =
      //       bmc.get<1>().find(boost::make_tuple(field, field, row_type, row_type));
      //   if (bit != bmc.get<1>().end()) {
 
      //     auto &m = *const_cast<MatrixDouble *>(&(bit->M));
      //     CHKERR invert_nonsymm_mat(m, inv);
 
      //     if (debug) {
      //       std::cerr << prod(m, inv) << endl;
      //       std::cerr << endl;
      //     }
 
      //   } else
      //     SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
      //              "%s matrix not found", field.c_str());
 
      //   MoFEMFunctionReturnHot(0);
      // };
 
      auto create_block_schur = [&](BlockMatContainer &bmc,
                                    BlockMatContainer &add_bmc,
                                    std::string field, AO ao,
                                    MatrixDouble &inv) {
        MoFEMFunctionBeginHot;
        // AOView(ao, PETSC_VIEWER_STDOUT_SELF);
        for (auto &bit : add_bmc) {
          bit.unSetAtElement();
          bit.clearMat();
        }
 
        for (auto &bit : bmc) {
          if (bit.setAtElement && bit.rowField != field &&
              bit.colField != field) {
            VectorInt rind = bit.rowInd;
            VectorInt cind = bit.colInd;
            const MatrixDouble &m = bit.M;
            if (ao) {
              CHKERR AOApplicationToPetsc(ao, rind.size(), &*rind.begin());
              CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
            }
            auto it =
                set_block(add_bmc, bit.rowField, bit.colField, bit.rowSide,
                          bit.colSide, bit.rowType, bit.colType, m, rind, cind);
          }
        }
 
        for (auto &bit_col : bmc) {
          if (bit_col.setAtElement && bit_col.rowField == field &&
              bit_col.colField != field) {
            const MatrixDouble &cm = bit_col.M;
            VectorInt cind = bit_col.colInd;
            invMat.resize(inv.size1(), cm.size2(), false);
            noalias(invMat) = prod(inv, cm);
            if (ao)
              CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
            for (auto &bit_row : bmc) {
              if (bit_row.setAtElement && bit_row.rowField != field &&
                  bit_row.colField == field) {
                const MatrixDouble &rm = bit_row.M;
                VectorInt rind = bit_row.rowInd;
                K.resize(rind.size(), cind.size(), false);
                noalias(K) = prod(rm, invMat);
                if (ao)
                  CHKERR AOApplicationToPetsc(ao, rind.size(), &*rind.begin());
                auto it =
                    add_block(add_bmc, bit_row.rowField, bit_col.colField,
                              bit_row.rowSide, bit_col.colSide, bit_row.rowType,
                              bit_col.colType, K, rind, cind);
              }
            }
          }
        }
 
        MoFEMFunctionReturnHot(0);
      };
 
      auto assemble_schur = [&](BlockMatContainer &add_bmc, Mat S,
                                bool debug = false) {
        MoFEMFunctionBeginHot;
        if (debug) {
          for (auto &bit : add_bmc) {
            if (bit.setAtElement) {
              std::cerr << "assemble vol: " << bit.rowField << " "
                        << bit.colField << endl;
              // std::cerr << bit.M << endl;
            }
          }
          std::cerr << std::endl;
        }
        for (auto &bit : add_bmc) {
          if (bit.setAtElement)
            CHKERR assemble_block(bit, S);
        }
        MoFEMFunctionReturnHot(0);
      };
 
      auto precondition_schur = [&](BlockMatContainer &bmc,
                                    BlockMatContainer &add_bmc,
                                    const std::string field,
                                    const MatrixDouble &diag_mat,
                                    const double eps) {
        MoFEMFunctionBeginHot;
 
        for (auto &bit : add_bmc) {
          bit.unSetAtElement();
          bit.clearMat();
        }
 
        for (auto &bit : bmc) {
          if (bit.setAtElement) {
            if (bit.rowField != field || bit.colField != field)
              auto it = set_block(add_bmc, bit.rowField, bit.colField,
                                  bit.rowSide, bit.colSide, bit.rowType,
                                  bit.colType, bit.M, bit.rowInd, bit.colInd);
          }
        }
 
        auto bit =
            bmc.get<1>().find(boost::make_tuple(field, field, row_type, row_type));
        if (bit->setAtElement && bit != bmc.get<1>().end()) {
          auto it = set_block(add_bmc, bit->rowField, bit->colField,
                              bit->rowSide, bit->colSide, bit->rowType,
                              bit->colType, bit->M, bit->rowInd, bit->colInd);
          MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
          m += eps * diag_mat;
        } else {
          auto row_it = bmc.get<3>().lower_bound(field);
          for (; row_it != bmc.get<3>().upper_bound(field); ++row_it) {
            if (row_it->setAtElement) {
              auto it = set_block(add_bmc, field, field, 0, 0, row_type, row_type,
                                  diag_mat, row_it->rowInd, row_it->rowInd);
              MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
              m *= eps;
              break;
            }
          }
          if (row_it == bmc.get<3>().end())
            SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
                     "row field not found %s", field.c_str());
        }
 
        MoFEMFunctionReturnHot(0);
      };
 
      const bool debug = false;
      auto &mat_cont = *commonDataPtr->blockMatContainerPtr;
 
      CHKERR invert_symm_schur(mat_cont, "EP",
                               invBlockMat["EPEP"]);
      CHKERR create_block_schur(mat_cont,
                                blockMat["EPEP"], "EP", aoSEpEp,
                                invBlockMat["EPEP"]);
      
      auto &mass_mat_scalar = commonDataPtr->massMatScalar;
      constexpr double eps = 1e-8; //FIXME: TODO: figure out this later
      CHKERR precondition_schur(blockMat["EPEP"], blockMat["precEPEP"], "TAU",
                                mass_mat_scalar, eps);
      // if (getPtrFE()->mField.check_field("SIGMA") && false) {
      //   auto &mass_mat_tensor = commonDataPtr->massMatTensor;
      //   CHKERR precondition_schur(blockMat["precEPEP"], blockMat["precSIGMA"],
      //                             "SIGMA", mass_mat_tensor, eps);
      //   CHKERR assemble_schur(blockMat["precSIGMA"], SEpEp, debug);
      //   blockMat["precEPEP"] = blockMat["precSIGMA"];
      // } else
      CHKERR assemble_schur(blockMat["precEPEP"], SEpEp, false);
 
      if(STauTau){
        CHKERR invert_symm_schur(blockMat["precEPEP"], "TAU",
                                 invBlockMat["TAUTAU"]);
        CHKERR create_block_schur(blockMat["precEPEP"], blockMat["TAUTAU"], "TAU",
                                  aoSTauTau, invBlockMat["TAUTAU"]);
        CHKERR assemble_schur(blockMat["TAUTAU"], STauTau);
      }
      
      blockMat.clear();
 
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
template struct OpInternalForceRhs<true, true>;
template struct OpInternalForceRhs<false, false>;
template struct OpInternalForceRhs<true, false>;
template struct OpInternalForceRhs<false, true>;
 
template struct OpPostProcElastic<true>;
template struct OpPostProcElastic<false>;
 
template struct OpCalculateLogStressTangent<true>;
template struct OpCalculateLogStressTangent<false>;
 
} // namespace OpElasticTools