fluid_structure_eigenproblem.cpp

Go to the documentation of this file.

/**
 * \file eigen_elastic.cpp
 * \FSI eigen_elastic.cpp
 *
 * Calculate natural frequencies in 2d and 3d problems.
 *
 */
 
#include <MoFEM.hpp>
#undef EPS
#include <slepceps.h>
#include <slepcpep.h>
 
using namespace MoFEM;
 
template <int DIM> struct ElementsAndOps {};
 
template <> struct ElementsAndOps<2> {
  using DomainEle = FaceElementForcesAndSourcesCore;
  using DomainEleOp = DomainEle::UserDataOperator;
  using PostProcEle = PostProcFaceOnRefinedMesh;
 
  using FaceSideEle = MoFEM::FaceElementForcesAndSourcesCoreOnSide;
  using FaceSideOp = FaceSideEle::UserDataOperator;
 
  using SkeletonEle = MoFEM::EdgeElementForcesAndSourcesCore;
  using SkeletonEleOp = SkeletonEle::UserDataOperator;
};
 
template <> struct ElementsAndOps<3> {
  using DomainEle = VolumeElementForcesAndSourcesCore;
  using DomainEleOp = DomainEle::UserDataOperator;
  using PostProcEle = PostProcVolumeOnRefinedMesh;
 
  using FaceSideEle = MoFEM::VolumeElementForcesAndSourcesCoreOnSide;
  using FaceSideOp = FaceSideEle::UserDataOperator;
 
  using SkeletonEle = MoFEM::FaceElementForcesAndSourcesCore;
  using SkeletonEleOp = SkeletonEle::UserDataOperator;
};
 
constexpr int SPACE_DIM =
    EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
 
using EntData = EntitiesFieldData::EntData;
using DomainEle = ElementsAndOps<SPACE_DIM>::DomainEle;
using DomainEleOp = ElementsAndOps<SPACE_DIM>::DomainEleOp;
using PostProcEle = ElementsAndOps<SPACE_DIM>::PostProcEle;
 
using FaceSideEle = ElementsAndOps<SPACE_DIM>::FaceSideEle;
using FaceSideOp = ElementsAndOps<SPACE_DIM>::FaceSideOp;
 
using SkeletonEle = ElementsAndOps<SPACE_DIM>::SkeletonEle;
using SkeletonEleOp = ElementsAndOps<SPACE_DIM>::SkeletonEleOp;
 
using OpK = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
    GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
using OpMass = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
    GAUSS>::OpMass<1, SPACE_DIM>;
 
FTensor::Index<'i', SPACE_DIM> i;
FTensor::Index<'j', SPACE_DIM> j;
FTensor::Index<'k', SPACE_DIM> k;
FTensor::Index<'l', SPACE_DIM> l;
FTensor::Index<'m', SPACE_DIM> m;
FTensor::Index<'n', SPACE_DIM> n;
 
static double penalty = 1e2;
static double phi = -1;
 
// Length unit are milometer
// Time unit miliseconds
// Force unit Newton
 
double rho_solid = 2.28;
double K_solid = 95 * 1e3;
double G_solid = 63 * 1e3;
 
double rho_fluid = 0.998;
double K_fluid = 2.2 * 1e3;
double G_fluid = 1e-6;
double MU_fluid = 1e-3;
 
int order = 1;
 
double eig_threshold = 1e-1;
 
// FIXME
auto poisson_ratio = 0.265;
#include <OpPostProcElastic.hpp>
 
using namespace Tutorial;
 
/**
 * @brief Operator tp collect data from elements on the side of Edge/Face
 *
 */
struct OpCalculateSideData : public FaceSideOp {
 
  OpCalculateSideData(std::string field_name, std::string col_field_name,
                      UserDataOperator::OpType type);
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
};
 
/**
 * @brief Operator the left hand side matrix
 *
 */
struct OpLhsSkeleton : public SkeletonEleOp {
public:
  /**
   * @brief Construct a new OpH1LhsSkeleton
   *
   * @param side_fe_ptr pointer to FE to evaluate side elements
   */
  OpLhsSkeleton(boost::shared_ptr<FaceSideEle> side_fe_ptr,
                boost::shared_ptr<MatrixDouble> d_mat_fluid_ptr,
                boost::shared_ptr<MatrixDouble> d_mat_solid_ptr,
                boost::shared_ptr<MatrixDouble> d_mat_visc_fluid_ptr,
                boost::shared_ptr<MatrixDouble> d_mat_visc_solid_ptr,
                boost::shared_ptr<Range> r_visc_fluid_ptr,
                boost::shared_ptr<Range> r_visc_solid_ptr, SmartPetscObj<Mat> K,
                SmartPetscObj<Mat> C, SmartPetscObj<Mat> M);
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data);
 
private:
  boost::shared_ptr<FaceSideEle>
      sideFEPtr; ///< pointer to element to get data on edge/face sides
  MatrixDouble locMatK;
  MatrixDouble locMatC;
  MatrixDouble locMatM;
  boost::shared_ptr<MatrixDouble> dMatFluidPtr;
  boost::shared_ptr<MatrixDouble> dMatSolidPtr;
  boost::shared_ptr<MatrixDouble> dViscFluidPtr;
  boost::shared_ptr<MatrixDouble> dViscSolidPtr;
  boost::shared_ptr<Range> rFluidPtr;
  boost::shared_ptr<Range> rSolidPtr;
  SmartPetscObj<Mat> K;
  SmartPetscObj<Mat> C;
  SmartPetscObj<Mat> M;
};
 
struct FSI {
 
  FSI(MoFEM::Interface &m_field) : mField(m_field) {}
 
  MoFEMErrorCode runProblem();
 
private:
  MoFEM::Interface &mField;
 
  MoFEMErrorCode readMesh();
  MoFEMErrorCode setupProblem();
  MoFEMErrorCode createCommonData();
  MoFEMErrorCode boundaryCondition();
  MoFEMErrorCode assembleSystem();
  MoFEMErrorCode solveSystem();
  MoFEMErrorCode outputResults();
 
  boost::shared_ptr<MatrixDouble> matGradPtr;
  boost::shared_ptr<MatrixDouble> matStrainPtr;
  boost::shared_ptr<MatrixDouble> matStressPtr;
  boost::shared_ptr<MatrixDouble> matDSolidPtr;
  boost::shared_ptr<MatrixDouble> matDFluidPtr;
  boost::shared_ptr<MatrixDouble> matDFViscFluidPtr;
  boost::shared_ptr<MatrixDouble> matDFViscSolidPtr;
 
  SmartPetscObj<Mat> M;
  SmartPetscObj<Mat> C;
  SmartPetscObj<Mat> K;
 
  SmartPetscObj<EPS> ePS;
  SmartPetscObj<PEP> pEP;
 
  enum SolverType { EPS_SOLVER, PEP_SOLVER, LASTSOLVER };
  int solverType;
 
  Range rSolid, rFluid, rInterface;
};
 
//! [Create common data]
MoFEMErrorCode FSI::createCommonData() {
  MoFEMFunctionBegin;
 
  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-rho_solid", &rho_solid,
                               PETSC_NULL);
  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-rho_fluid", &rho_fluid,
                               PETSC_NULL);
  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-bulk_modulus_solid", &K_solid,
                               PETSC_NULL);
  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-sheer_modulus_solid", &G_solid,
                               PETSC_NULL);
  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-bulk_modulus_fluid", &K_fluid,
                               PETSC_NULL);
  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-viscosity_fluid", &MU_fluid,
                               PETSC_NULL);
  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-penalty", &penalty,
                               PETSC_NULL);
  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-phi", &phi, PETSC_NULL);
 
  MOFEM_LOG("FSI", Sev::inform) << "-rho_solid " << rho_solid;
  MOFEM_LOG("FSI", Sev::inform) << "-rho_fluid " << rho_fluid;
  MOFEM_LOG("FSI", Sev::inform) << "-bulk_modulus_solid " << K_solid;
  MOFEM_LOG("FSI", Sev::inform) << "-sheer_modulus_solid " << G_solid;
  MOFEM_LOG("FSI", Sev::inform) << "-bulk_modulus_fluid " << K_fluid;
  MOFEM_LOG("FSI", Sev::inform) << "-viscosity_fluid " << MU_fluid;
  MOFEM_LOG("FSI", Sev::inform) << "-penalty " << penalty;
  MOFEM_LOG("FSI", Sev::inform) << "-phi " << phi;
 
  const char *list_solvers[]{"eps", "pep"};
  solverType = EPS_SOLVER;
  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-solver_type", list_solvers,
                              LASTSOLVER, &solverType, PETSC_NULL);
 
  auto set_matrial_stiffens = [&](auto &mat_d, auto K, auto G) {
    constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
    mat_d.resize(size_symm * size_symm, 1);
 
    constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
    MoFEMFunctionBegin;
    auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(mat_d);
 
    constexpr double A = 1; // plane strain
 
    // double A = (SPACE_DIM == 2) ? 2 * G / (K + (4. / 3.) * G) : 1;
 
    t_D(i, j, k, l) = 2 * G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
                      A * (K - (2. / 3.) * G) * t_kd(i, j) * t_kd(k, l);
 
    MoFEMFunctionReturn(0);
  };
 
  auto set_viscosity = [&](auto &mat_d, auto MU) {
    constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
    mat_d.resize(size_symm * size_symm, 1);
 
    constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
    MoFEMFunctionBegin;
    auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(mat_d);
 
    t_D(i, j, k, l) = MU * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
 
    MoFEMFunctionReturn(0);
  };
 
  matGradPtr = boost::make_shared<MatrixDouble>();
  matStrainPtr = boost::make_shared<MatrixDouble>();
  matStressPtr = boost::make_shared<MatrixDouble>();
  matDSolidPtr = boost::make_shared<MatrixDouble>();
  matDFluidPtr = boost::make_shared<MatrixDouble>();
  matDFViscFluidPtr = boost::make_shared<MatrixDouble>();
  matDFViscSolidPtr = boost::make_shared<MatrixDouble>();
 
  CHKERR set_matrial_stiffens(*matDSolidPtr, K_solid, G_solid);
  CHKERR set_matrial_stiffens(*matDFluidPtr, K_fluid, G_fluid);
  // CHKERR set_matrial_stiffens(*matDFViscFluidPtr, 0, MU_fluid);
  CHKERR set_viscosity(*matDFViscFluidPtr, MU_fluid);
 
  // No solid damping
  matDFViscSolidPtr->resize(matDFViscFluidPtr->size1(),
                            matDFViscFluidPtr->size2());
  matDFViscSolidPtr->clear();
 
  MoFEMFunctionReturn(0);
}
//! [Create common data]
 
//! [Run problem]
MoFEMErrorCode FSI::runProblem() {
  MoFEMFunctionBegin;
  CHKERR createCommonData();
  CHKERR readMesh();
  CHKERR setupProblem();
  CHKERR boundaryCondition();
  CHKERR assembleSystem();
  CHKERR solveSystem();
  CHKERR outputResults();
  MoFEMFunctionReturn(0);
}
//! [Run problem]
 
//! [Read mesh]
MoFEMErrorCode FSI::readMesh() {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
 
  MOFEM_LOG("FSI", Sev::inform)
      << "Read mesh for problem in " << EXECUTABLE_DIMENSION;
  CHKERR simple->getOptions();
 
  simple->getAddSkeletonFE() = true;
  CHKERR simple->loadFile();
 
  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
 
    const std::string block_names[] = {"SOLID", "FLUID", "INTERFACE"};
 
    if (it->getName().compare(0, block_names[0].length(), block_names[0]) ==
        0) {
      CHKERR mField.get_moab().get_entities_by_handle(it->meshset, rSolid,
                                                      true);
      rSolid = rSolid.subset_by_dimension(SPACE_DIM);
    } else if (it->getName().compare(0, block_names[1].length(),
                                     block_names[1]) == 0) {
      CHKERR mField.get_moab().get_entities_by_handle(it->meshset, rFluid,
                                                      true);
      rFluid = rFluid.subset_by_dimension(SPACE_DIM);
    } else if (it->getName().compare(0, block_names[2].length(),
                                     block_names[2]) == 0) {
      CHKERR mField.get_moab().get_entities_by_handle(it->meshset, rInterface,
                                                      true);
      rInterface = rInterface.subset_by_dimension(SPACE_DIM - 1);
    }
  }
 
  for (auto d = SPACE_DIM - 1; d >= 0; --d) {
    auto add_adj = [&](auto &r, auto dim) {
      MoFEMFunctionBeginHot;
      Range adj;
      CHKERR mField.get_moab().get_adjacencies(
          r.subset_by_dimension(dim), d, false, adj, moab::Interface::UNION);
      r.merge(adj);
      MoFEMFunctionReturnHot(0);
    };
    CHKERR add_adj(rSolid, SPACE_DIM);
    CHKERR add_adj(rFluid, SPACE_DIM);
    if (d < (SPACE_DIM - 1)) {
      CHKERR add_adj(rInterface, SPACE_DIM - 1);
    }
  }
 
  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(rSolid);
  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(rFluid);
  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(rInterface);
 
  MOFEM_LOG("FSI", Sev::noisy) << "Solid " << rSolid << endl;
  MOFEM_LOG("FSI", Sev::noisy) << "Fluid " << rFluid << endl;
  MOFEM_LOG("FSI", Sev::noisy) << "Interface " << rInterface << endl;
 
  MoFEMFunctionReturn(0);
}
//! [Read mesh]
 
//! [Set up problem]
MoFEMErrorCode FSI::setupProblem() {
  auto *simple = mField.getInterface<Simple>();
  MoFEMFunctionBegin;
  // Add field
  CHKERR simple->addDomainField("Us", H1, AINSWORTH_LEGENDRE_BASE, SPACE_DIM);
  CHKERR simple->addDomainField("Uf", H1, AINSWORTH_LEGENDRE_BASE, SPACE_DIM);
  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
  CHKERR simple->setFieldOrder("Us", order);
  CHKERR simple->setFieldOrder("Uf", order);
 
  if (solverType == PEP_SOLVER) {
 
    CHKERR mField.add_finite_element("INTERFACE");
    auto add_fe_field = [&](auto field) {
      MoFEMFunctionBeginHot;
      CHKERR mField.modify_finite_element_add_field_row("INTERFACE", field);
      CHKERR mField.modify_finite_element_add_field_col("INTERFACE", field);
      CHKERR mField.modify_finite_element_add_field_data("INTERFACE", field);
      MoFEMFunctionReturnHot(0);
    };
    CHKERR add_fe_field("Uf");
    CHKERR add_fe_field("Us");
 
    CHKERR mField.add_ents_to_finite_element_by_dim(rInterface, SPACE_DIM - 1,
                                                    "INTERFACE");
    simple->getOtherFiniteElements().push_back("INTERFACE");
  }
 
  CHKERR simple->defineFiniteElements();
  CHKERR simple->setSkeletonAdjacency();
  CHKERR simple->defineProblem();
  CHKERR simple->buildFields();
  CHKERR simple->buildFiniteElements();
  CHKERR simple->buildProblem();
 
  CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
      simple->getProblemName(), "Uf", subtract(rSolid, rInterface));
  CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
      simple->getProblemName(), "Us", subtract(rFluid, rInterface));
 
  MoFEMFunctionReturn(0);
}
//! [Set up problem]
 
//! [Boundary condition]
MoFEMErrorCode FSI::boundaryCondition() {
  auto simple = mField.getInterface<Simple>();
  auto bc_mng = mField.getInterface<BcManager>();
  MoFEMFunctionBegin;
 
  auto bc_mng = mField.getInterface<BcManager>();
  auto &bc_map = bc_mng->getBcMapByBlockName();
 
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_X",
                                           "Uf", 0, 0);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Y",
                                           "Uf", 1, 1);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Z",
                                           "Uf", 2, 2);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
                                           "Uf", 0, 3);
 
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_X",
                                           "Us", 0, 0);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Y",
                                           "Us", 1, 1);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Z",
                                           "Us", 2, 2);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
                                           "Us", 0, 3);
 
  MoFEMFunctionReturn(0);
}
//! [Boundary condition]
 
//! [Push operators to pipeline]
MoFEMErrorCode FSI::assembleSystem() {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto pipeline_mng = mField.getInterface<PipelineManager>();
  auto bc_mng = mField.getInterface<BcManager>();
  auto &bc_map = bc_mng->getBcMapByBlockName();
 
  auto dm = simple->getDM();
  CHKERR DMCreateMatrix_MoFEM(dm, K);
  M = smartMatDuplicate(K, MAT_SHARE_NONZERO_PATTERN);
  C = smartMatDuplicate(K, MAT_SHARE_NONZERO_PATTERN);
 
  CHKERR MatZeroEntries(K);
  CHKERR MatZeroEntries(M);
  CHKERR MatZeroEntries(C);
 
  auto ptr_r_solid =
      boost::make_shared<Range>(rSolid.subset_by_dimension(SPACE_DIM));
  auto ptr_r_fluid =
      boost::make_shared<Range>(rFluid.subset_by_dimension(SPACE_DIM));
 
  auto basic_ops = [&]() {
    MoFEMFunctionBegin;
    auto det_ptr = boost::make_shared<VectorDouble>();
    auto jac_ptr = boost::make_shared<MatrixDouble>();
    auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
    pipeline_mng->getOpDomainLhsPipeline().push_back(
        new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
    pipeline_mng->getOpDomainLhsPipeline().push_back(
        new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
    pipeline_mng->getOpDomainLhsPipeline().push_back(
        new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
    pipeline_mng->getOpDomainLhsPipeline().push_back(
        new OpSetHOWeightsOnFace());
    MoFEMFunctionReturn(0);
  };
 
  auto integration_rule = [](int, int, int approx_order) {
    return 2 * approx_order + 1;
  };
 
  auto calculate_stiffness_matrix = [&]() {
    MoFEMFunctionBegin;
    pipeline_mng->getDomainLhsFE().reset();
    CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
 
    CHKERR basic_ops();
 
    pipeline_mng->getOpDomainLhsPipeline().push_back(
        new OpK("Us", "Us", matDSolidPtr, ptr_r_solid));
    pipeline_mng->getOpDomainLhsPipeline().push_back(
        new OpK("Uf", "Uf", matDFluidPtr, ptr_r_fluid));
 
    pipeline_mng->getDomainLhsFE()->B = K;
    CHKERR pipeline_mng->loopFiniteElements();
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_L_matrix = [&]() {
    MoFEMFunctionBegin;
 
    auto get_side_fe = [&]() {
      auto side_fe_ptr = boost::make_shared<FaceSideEle>(mField);
      auto det_ptr = boost::make_shared<VectorDouble>();
      auto jac_ptr = boost::make_shared<MatrixDouble>();
      auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
      side_fe_ptr->getOpPtrVector().push_back(
          new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
      side_fe_ptr->getOpPtrVector().push_back(
          new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
      side_fe_ptr->getOpPtrVector().push_back(
          new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
      side_fe_ptr->getOpPtrVector().push_back(
          new OpCalculateSideData("Us", "Us", FaceSideOp::OPROW));
      side_fe_ptr->getOpPtrVector().push_back(
          new OpCalculateSideData("Uf", "Uf", FaceSideOp::OPCOL));
      return side_fe_ptr;
    };
 
    auto fe_side_ptr = boost::make_shared<SkeletonEle>(mField);
    fe_side_ptr->getOpPtrVector().push_back(new OpLhsSkeleton(
        get_side_fe(), matDFluidPtr, matDSolidPtr, matDFViscFluidPtr,
        matDFViscSolidPtr, boost::make_shared<Range>(rFluid),
        boost::make_shared<Range>(rSolid), K, C, M));
 
    auto rules = fe_side_ptr->getRuleHook = [](int, int, int p) -> int {
      return 2 * p + 1;
    };
 
    CHKERR mField.loop_finite_elements(simple->getProblemName(), "INTERFACE",
                                       *fe_side_ptr);
 
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_damping_matrix = [&]() {
    MoFEMFunctionBegin;
 
    pipeline_mng->getDomainLhsFE().reset();
    CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
 
    CHKERR basic_ops();
 
    pipeline_mng->getOpDomainLhsPipeline().push_back(
        new OpK("Uf", "Uf", matDFViscFluidPtr, ptr_r_fluid));
 
    pipeline_mng->getDomainLhsFE()->B = C;
    CHKERR pipeline_mng->loopFiniteElements();
 
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_mass_matrix = [&]() {
    MoFEMFunctionBegin;
 
    pipeline_mng->getDomainLhsFE().reset();
    CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
 
    CHKERR basic_ops();
 
    pipeline_mng->getOpDomainLhsPipeline().push_back(new OpMass(
        "Us", "Us",
        [](const double, const double, const double) { return rho_solid; },
        ptr_r_solid));
    pipeline_mng->getOpDomainLhsPipeline().push_back(new OpMass(
        "Uf", "Uf",
        [](const double, const double, const double) { return rho_fluid; },
        ptr_r_fluid));
 
    pipeline_mng->getDomainLhsFE()->B = M;
    CHKERR pipeline_mng->loopFiniteElements();
    MoFEMFunctionReturn(0);
  };
 
  if (solverType == PEP_SOLVER)
    CHKERR calculate_L_matrix();
 
  CHKERR calculate_stiffness_matrix();
  CHKERR calculate_damping_matrix();
  CHKERR calculate_mass_matrix();
 
  CHKERR MatAssemblyBegin(K, MAT_FINAL_ASSEMBLY);
  CHKERR MatAssemblyEnd(K, MAT_FINAL_ASSEMBLY);
  CHKERR MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
  CHKERR MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
  CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
  CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
 
  // auto save_range = [&](const std::string name, const Range r) {
  //   MoFEMFunctionBegin;
  //   EntityHandle out_meshset;
  //   CHKERR mField.get_moab().create_meshset(MESHSET_SET, out_meshset);
  //   CHKERR mField.get_moab().add_entities(out_meshset, r);
  //   CHKERR mField.get_moab().write_file(name.c_str(), "VTK", "",
  //   &out_meshset,
  //                                       1);
  //   CHKERR mField.get_moab().delete_entities(&out_meshset, 1);
  //   MoFEMFunctionReturn(0);
  // };
 
  // CHKERR save_range("interface.vtk", rInterface);
 
  MoFEMFunctionReturn(0);
}
//! [Push operators to pipeline]
 
//! [Solve]
MoFEMErrorCode FSI::solveSystem() {
  auto simple = mField.getInterface<Simple>();
  auto is_manager = mField.getInterface<ISManager>();
  MoFEMFunctionBegin;
 
  auto create_eps = [&](MPI_Comm comm) {
    EPS eps;
    CHKERR EPSCreate(comm, &eps);
    CHKERR EPSSetOperators(eps, K, M);
    return SmartPetscObj<EPS>(eps);
  };
 
  auto print_info_eps = [&](auto eps) {
    MoFEMFunctionBegin;
    ST st;
    EPSType type;
    PetscReal tol;
    PetscInt nev, maxit, its;
    // Optional: Get some information from the solver and display it
    CHKERR EPSGetIterationNumber(eps, &its);
    MOFEM_LOG_C("FSI", Sev::inform, " Number of iterations of the method: %d",
                its);
    CHKERR EPSGetST(eps, &st);
    CHKERR EPSGetType(eps, &type);
    MOFEM_LOG_C("FSI", Sev::inform, " Solution method: %s", type);
    CHKERR EPSGetDimensions(eps, &nev, NULL, NULL);
    MOFEM_LOG_C("FSI", Sev::inform, " Number of requested eigenvalues: %d",
                nev);
    CHKERR EPSGetTolerances(eps, &tol, &maxit);
    MOFEM_LOG_C("FSI", Sev::inform, " Stopping condition: tol=%.4g, maxit=%d",
                (double)tol, maxit);
 
    PetscScalar eigr, eigi;
    for (int nn = 0; nn < nev; nn++) {
      CHKERR EPSGetEigenpair(eps, nn, &eigr, &eigi, PETSC_NULL, PETSC_NULL);
      MOFEM_LOG_C("FSI", Sev::inform,
                  " ncov = %d eigr = %.8g eigi = %.8g (inv eigr = %.8g)", nn,
                  eigr, eigi, 1. / eigr);
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto setup_eps = [&](auto eps) {
    MoFEMFunctionBegin;
    CHKERR EPSSetFromOptions(eps);
    MoFEMFunctionReturn(0);
  };
 
  Mat A[3] = {K.get(), C.get(), M.get()};
 
  auto create_pep = [&](MPI_Comm comm) {
    PEP pep;
    CHKERR PEPCreate(comm, &pep);
    CHKERR PEPSetOperators(pep, 3, A);
    return SmartPetscObj<PEP>(pep);
  };
 
  auto setup_pep = [&](auto pep) {
    MoFEMFunctionBegin;
    CHKERR PEPSetFromOptions(pep);
    MoFEMFunctionReturn(0);
  };
 
  auto print_info_pep = [&](auto pep) {
    MoFEMFunctionBegin;
    int nconv;
    double kr, ki, error;
    CHKERR PEPGetConverged(pEP, &nconv);
    for (int j = 0; j < nconv; ++j) {
      CHKERR PEPGetEigenpair(pep, j, &kr, &ki, PETSC_NULL, PETSC_NULL);
      MOFEM_LOG_C("FSI", Sev::verbose, " ncov = %d eigr = %.8g eigi = %.8g", j,
                  kr, ki);
    }
 
    MoFEMFunctionReturn(0);
  };
 
  if (solverType == EPS_SOLVER) {
    // Create eigensolver context
    ePS = create_eps(mField.get_comm());
    // Setup eps
    CHKERR setup_eps(ePS);
    // Solve problem
    CHKERR EPSSolve(ePS);
    // Print info
    CHKERR print_info_eps(ePS);
  } else if (solverType == PEP_SOLVER) {
    pEP = create_pep(mField.get_comm());
    CHKERR setup_pep(pEP);
    CHKERR PEPSolve(pEP);
    CHKERR print_info_pep(pEP);
  } else {
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Solver not implemented");
  }
 
  MoFEMFunctionReturn(0);
}
//! [Solve]
 
//! [Postprocess results]
MoFEMErrorCode FSI::outputResults() {
  MoFEMFunctionBegin;
  auto *pipeline_mng = mField.getInterface<PipelineManager>();
  auto *simple = mField.getInterface<Simple>();
 
  pipeline_mng->getDomainLhsFE().reset();
  pipeline_mng->getDomainRhsFE().reset();
  auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
  post_proc_fe->generateReferenceElementMesh();
 
  auto det_ptr = boost::make_shared<VectorDouble>();
  auto jac_ptr = boost::make_shared<MatrixDouble>();
  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
  post_proc_fe->getOpPtrVector().push_back(
      new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
  post_proc_fe->getOpPtrVector().push_back(
      new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
  post_proc_fe->getOpPtrVector().push_back(
      new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
 
  auto us_ptr = boost::make_shared<MatrixDouble>();
  auto grad_ptr = boost::make_shared<MatrixDouble>();
  auto strain_s_ptr = boost::make_shared<MatrixDouble>();
  auto uf_ptr = boost::make_shared<MatrixDouble>();
  auto strain_f_ptr = boost::make_shared<MatrixDouble>();
 
  // auto stress_ptr = boost::make_shared<MatrixDouble>();
 
  post_proc_fe->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<SPACE_DIM>("Us", us_ptr));
  post_proc_fe->getOpPtrVector().push_back(
      new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>("Us", grad_ptr));
  post_proc_fe->getOpPtrVector().push_back(
      new OpSymmetrizeTensor<SPACE_DIM>("Us", grad_ptr, strain_s_ptr));
 
  post_proc_fe->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<SPACE_DIM>("Uf", uf_ptr));
  post_proc_fe->getOpPtrVector().push_back(
      new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>("Uf", grad_ptr));
  post_proc_fe->getOpPtrVector().push_back(
      new OpSymmetrizeTensor<SPACE_DIM>("Uf", grad_ptr, strain_f_ptr));
 
  using OpPPMap = OpPostProcMap<SPACE_DIM, SPACE_DIM>;
 
  post_proc_fe->getOpPtrVector().push_back(
 
      new OpPPMap(post_proc_fe->postProcMesh, post_proc_fe->mapGaussPts,
 
                  OpPPMap::DataMapVec{},
 
                  OpPPMap::DataMapMat{{"Us", us_ptr}, {"Uf", uf_ptr}},
 
                  OpPPMap::DataMapMat{},
 
                  OpPPMap::DataMapMat{{"STRAINs", strain_s_ptr},
                                      {"STRAINf", strain_f_ptr}}
 
                  )
 
  );
 
  auto bc_mng = mField.getInterface<BcManager>();
  auto &bc_map = bc_mng->getBcMapByBlockName();
 
  auto op_tag_handle = new DomainEleOp(NOSPACE, DomainEleOp::OPSPACE);
  op_tag_handle->doWorkRhsHook = [&](DataOperator *op_ptr, int side,
                                     EntityType type,
                                     EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
 
    auto get_tag = [&]() {
      int def_val = -1;
      Tag th;
      CHK_MOAB_THROW(post_proc_fe->postProcMesh.tag_get_handle(
                         "BLOCK", 1, MB_TYPE_INTEGER, th,
                         MB_TAG_CREAT | MB_TAG_SPARSE, &def_val),
                     "Can not create tag on post-proc mesh");
      return th;
    };
 
    auto fe_ent = static_cast<DomainEleOp *>(op_ptr)->getFEEntityHandle();
 
    int marker = 0;
    if (rFluid.find(fe_ent) != rFluid.end()) {
      marker += 1;
    }
    if (rSolid.find(fe_ent) != rSolid.end()) {
      marker += 2;
    }
 
    auto &nodes_vec = post_proc_fe->mapGaussPts;
    CHK_MOAB_THROW(
        post_proc_fe->postProcMesh.tag_clear_data(
            get_tag(), &*nodes_vec.begin(), nodes_vec.size(), &marker),
        "Can not set tag");
 
    MoFEMFunctionReturn(0);
  };
 
  post_proc_fe->getOpPtrVector().push_back(op_tag_handle);
 
  pipeline_mng->getDomainRhsFE() = post_proc_fe;
 
  auto dm = simple->getDM();
 
  auto save_vec = [&](auto file_name, auto v) {
    MoFEMFunctionBegin;
    CHKERR DMoFEMMeshToLocalVector(dm, v, INSERT_VALUES, SCATTER_REVERSE);
    CHKERR pipeline_mng->loopFiniteElements();
    post_proc_fe->writeFile(file_name);
    MoFEMFunctionReturn(0);
  };
 
  auto post_proc_eps = [&](auto eps) {
    MoFEMFunctionBegin;
    auto v = smartCreateDMVector(dm);
 
    PetscInt nev;
    CHKERR EPSGetDimensions(ePS, &nev, NULL, NULL);
    PetscScalar eigr, eigi, nrm2r;
    for (int nn = 0; nn < nev; nn++) {
      CHKERR EPSGetEigenpair(ePS, nn, &eigr, &eigi, v, PETSC_NULL);
      CHKERR VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecNorm(v, NORM_2, &nrm2r);
      MOFEM_LOG_C(
          "FSI", Sev::verbose,
          " ncov = %d omega2 = %.8g omega = %.8g osc frequency  = %.8g kHz",
          nn, eigr, std::sqrt(std::abs(eigr)),
          std::sqrt(std::abs(eigr)) / (2 * M_PI));
 
      CHKERR save_vec(
          "out_eig_" + boost::lexical_cast<std::string>(nn) + ".h5m", v);
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto post_proc_pep = [&](auto eps) {
    MoFEMFunctionBegin;
 
    auto re_vec = smartCreateDMVector(dm);
    auto im_vec = smartVectorDuplicate(re_vec);
 
    PetscInt nev;
    CHKERR PEPGetConverged(pEP, &nev);
    PetscScalar ki, kr, error, nrm2r;
 
    int idx = 0; // counting from zero, increment first before log
    for (int nn = 0; nn < nev; nn++) {
 
      CHKERR PEPGetEigenpair(pEP, nn, &kr, &ki, re_vec, im_vec);
      CHKERR PEPComputeError(pEP, nn, PEP_ERROR_BACKWARD, &error);
 
      CHKERR VecGhostUpdateBegin(re_vec, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(re_vec, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateBegin(im_vec, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(im_vec, INSERT_VALUES, SCATTER_FORWARD);
 
      if (std::abs(ki) > eig_threshold &&
          kr < std::numeric_limits<float>::epsilon()) {
        CHKERR save_vec("out_pep_re_" + boost::lexical_cast<std::string>(idx) +
                            ".h5m",
                        re_vec);
        CHKERR save_vec("out_pep_im_" + boost::lexical_cast<std::string>(idx) +
                            ".h5m",
                        im_vec);
        ++idx;
      }
 
      MOFEM_LOG_C("FSI", Sev::inform,
                  "idx = %d ncov = %d eigr = %.8g eigi = %.8g angular "
                  "frequency = %.8g kHz "
                  "error %.8g",
                  idx - 1, nn, kr, ki, ki / (2 * M_PI), error);
    }
    MoFEMFunctionReturn(0);
  };
 
  if (solverType == EPS_SOLVER)
    CHKERR post_proc_eps(ePS);
  else if (solverType == PEP_SOLVER)
    CHKERR post_proc_pep(pEP);
 
  MoFEMFunctionReturn(0);
}
//! [Postprocess results]
 
static char help[] = "...\n\n";
 
int main(int argc, char *argv[]) {
 
  // Initialisation of MoFEM/PETSc and MOAB data structures
  const char param_file[] = "param_file.petsc";
  SlepcInitialize(&argc, &argv, param_file, help);
  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
 
  // Add logging channel for FSI
  auto core_log = logging::core::get();
  core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "FSI"));
  LogManager::setLog("FSI");
  MOFEM_LOG_TAG("FSI", "FSI");
 
  try {
 
    //! [Register MoFEM discrete manager in PETSc]
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
    //! [Register MoFEM discrete manager in PETSc
 
    //! [Create MoAB]
    moab::Core mb_instance;              ///< mesh database
    moab::Interface &moab = mb_instance; ///< mesh database interface
    //! [Create MoAB]
 
    //! [Create MoFEM]
    MoFEM::Core core(moab);           ///< finite element database
    MoFEM::Interface &m_field = core; ///< finite element database insterface
    //! [Create MoFEM]
 
    //! [FSI]
    FSI ex(m_field);
    CHKERR ex.runProblem();
    //! [FSI]
  }
  CATCH_ERRORS;
 
  SlepcFinalize();
  CHKERR MoFEM::Core::Finalize();
}
 
// data for skeleton computation
std::array<std::vector<VectorInt>, 2> indicesUsSideMap;
std::array<std::vector<VectorInt>, 2> indicesUfSideMap;
 
std::array<std::vector<MatrixDouble>, 2> baseUfSideMap;
std::array<std::vector<MatrixDouble>, 2> diffUfBaseSideMap;
std::array<std::vector<MatrixDouble>, 2> baseUsSideMap;
std::array<std::vector<MatrixDouble>, 2> diffUsBaseSideMap;
 
std::array<EntityHandle, 2> sideHandle;
std::array<double, 2> areaMap; // area/volume of elements on the side
std::array<int, 2> senseMap;   // orientation of local element edge/face in
                               // respect to global orientation of edge/face
 
OpCalculateSideData::OpCalculateSideData(std::string row_field_name,
                                         std::string col_field_name,
                                         UserDataOperator::OpType type)
    : FaceSideOp(row_field_name, col_field_name, type) {}
 
MoFEMErrorCode OpCalculateSideData::doWork(int side, EntityType type,
                                           EntData &data) {
  MoFEMFunctionBeginHot;
 
  const auto nb_in_loop = getFEMethod()->nInTheLoop;
 
  auto clear = [](auto nb) {
    indicesUsSideMap[nb].clear();
    indicesUfSideMap[nb].clear();
    baseUsSideMap[nb].clear();
    diffUsBaseSideMap[nb].clear();
    baseUfSideMap[nb].clear();
    diffUfBaseSideMap[nb].clear();
  };
 
  if (getOpType() == OPROW) {
 
    if (type == MBVERTEX) {
      areaMap[nb_in_loop] = getMeasure();
      senseMap[nb_in_loop] = getSkeletonSense();
      sideHandle[nb_in_loop] = getFEEntityHandle();
      if (!nb_in_loop) {
        clear(0);
        clear(1);
        sideHandle[1] = 0;
        areaMap[1] = 0;
        senseMap[1] = 0;
      }
    }
 
    const auto nb_dofs = data.getIndices().size();
    if (nb_dofs) {
      indicesUsSideMap[nb_in_loop].push_back(data.getIndices());
      baseUsSideMap[nb_in_loop].push_back(
          data.getN(BaseDerivatives::ZeroDerivative));
      diffUsBaseSideMap[nb_in_loop].push_back(
          data.getN(BaseDerivatives::FirstDerivative));
    }
  } else if (getOpType() == OPCOL) {
 
    if (type == MBVERTEX) {
      if (areaMap[nb_in_loop] != getMeasure()) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Should be the same");
      }
      if (senseMap[nb_in_loop] != getSkeletonSense()) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Should be the same");
      }
      if (sideHandle[nb_in_loop] != getFEEntityHandle()) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Should be the same");
      };
    }
 
    const auto nb_dofs = data.getIndices().size();
    if (nb_dofs) {
      indicesUfSideMap[nb_in_loop].push_back(data.getIndices());
      baseUfSideMap[nb_in_loop].push_back(
          data.getN(BaseDerivatives::ZeroDerivative));
      diffUfBaseSideMap[nb_in_loop].push_back(
          data.getN(BaseDerivatives::FirstDerivative));
    }
 
  } else {
    SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong OpType %s",
             OpTypeNames[getOpType()]);
  }
 
  MoFEMFunctionReturnHot(0);
}
 
enum ElementSide { LEFT_SIDE = 0, RIGHT_SIDE };
 
template <typename T> inline auto get_ntensor(T &base_mat) {
  return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(
      &*base_mat.data().begin());
};
 
template <typename T> inline auto get_ntensor(T &base_mat, int gg, int bb) {
  double *ptr = &base_mat(gg, bb);
  return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(ptr);
};
 
template <typename T> inline auto get_diff_ntensor(T &base_mat) {
  double *ptr = &*base_mat.data().begin();
  return getFTensor1FromPtr<SPACE_DIM>(ptr);
};
 
template <typename T>
inline auto get_diff_ntensor(T &base_mat, int gg, int bb) {
  double *ptr = &base_mat(gg, SPACE_DIM * bb);
  return getFTensor1FromPtr<SPACE_DIM>(ptr);
};
 
/**
 * @brief Construct a new OpH1LhsSkeleton
 *
 * @param side_fe_ptr pointer to FE to evaluate side elements
 */
OpLhsSkeleton::OpLhsSkeleton(
    boost::shared_ptr<FaceSideEle> side_fe_ptr,
    boost::shared_ptr<MatrixDouble> mat_d_fluid_ptr,
    boost::shared_ptr<MatrixDouble> mat_d_solid_ptr,
    boost::shared_ptr<MatrixDouble> mat_d_visc_fluid_ptr,
    boost::shared_ptr<MatrixDouble> mat_d_visc_solid_ptr,
    boost::shared_ptr<Range> r_fluid_ptr, boost::shared_ptr<Range> r_solid_ptr,
    SmartPetscObj<Mat> K, SmartPetscObj<Mat> C, SmartPetscObj<Mat> M)
    : SkeletonEleOp(NOSPACE, SkeletonEleOp::OPSPACE), sideFEPtr(side_fe_ptr),
      dMatFluidPtr(mat_d_fluid_ptr), dMatSolidPtr(mat_d_solid_ptr),
      dViscFluidPtr(mat_d_visc_fluid_ptr), dViscSolidPtr(mat_d_visc_solid_ptr),
      rFluidPtr(r_fluid_ptr), rSolidPtr(r_solid_ptr), K(K), C(C), M(M) {}
 
MoFEMErrorCode OpLhsSkeleton::doWork(int side, EntityType type,
                                     EntitiesFieldData::EntData &data) {
  MoFEMFunctionBegin;
 
  // get normal of the face or edge
  auto t_normal = getFTensor1Normal();
  t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
 
  // Collect data from side domian elements
  CHKERR loopSide("dFE", sideFEPtr.get(), SPACE_DIM);
  const auto in_the_loop =
      sideFEPtr->nInTheLoop; // return number of elements on the side
 
  std::array<boost::shared_ptr<MatrixDouble>, 2> d_ptr;
  std::array<boost::shared_ptr<MatrixDouble>, 2> d_vis_ptr;
 
  std::array<std::vector<VectorInt>, 2> indices_side_map;
  std::array<std::vector<MatrixDouble>, 2> base_side_map;
  std::array<std::vector<MatrixDouble>, 2> diff_base_side_map;
 
  // MOFEM_LOG("SELF", Sev::inform)
  //     << "AAA " << sideHandle[0] << " " << sideHandle[1];
 
  // auto save_range = [&](const std::string name, const Range r) {
  //   MoFEMFunctionBegin;
  //   EntityHandle out_meshset;
  //   CHKERR getPtrFE()->mField.get_moab().create_meshset(MESHSET_SET,
  //                                                       out_meshset);
  //   CHKERR getPtrFE()->mField.get_moab().add_entities(out_meshset, r);
  //   CHKERR getPtrFE()->mField.get_moab().write_file(name.c_str(), "VTK", "",
  //                                                   &out_meshset, 1);
  //   CHKERR getPtrFE()->mField.get_moab().delete_entities(&out_meshset, 1);
  //   MoFEMFunctionReturn(0);
  // };
 
  // Tag th_n;
  // double def[] = {0, 0, 0};
  // CHKERR getPtrFE()->mField.get_moab().tag_get_handle(
  //     "NORMAL", 3, MB_TYPE_DOUBLE, th_n, MB_TAG_CREAT | MB_TAG_SPARSE, def);
 
  auto set_row_col = [&](auto a, auto b) {
    d_ptr[a] = dMatSolidPtr;
    d_ptr[b] = dMatFluidPtr;
    d_vis_ptr[a] = dViscSolidPtr;
    d_vis_ptr[b] = dViscFluidPtr;
    indices_side_map[a].swap(indicesUsSideMap[a]);
    indices_side_map[b].swap(indicesUfSideMap[b]);
    base_side_map[a].swap(baseUsSideMap[a]);
    base_side_map[b].swap(baseUfSideMap[b]);
    diff_base_side_map[a].swap(diffUsBaseSideMap[a]);
    diff_base_side_map[b].swap(diffUfBaseSideMap[b]);
  };
 
  if (rSolidPtr->find(sideHandle[0]) != rSolidPtr->end()) {
 
    if (rFluidPtr->find(sideHandle[1]) == rFluidPtr->end()) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Other side should be fluid");
    }
 
    set_row_col(0, 1);
 
    if (senseMap[0] < 0)
      t_normal(i) *= -1;
 
  } else if (rFluidPtr->find(sideHandle[0]) != rFluidPtr->end()) {
 
    if (rSolidPtr->find(sideHandle[1]) == rSolidPtr->end()) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Other side should be solid");
    }
 
    set_row_col(1, 0);
 
    if (senseMap[1] < 0)
      t_normal(i) *= -1;
 
  } else {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Side is not fluid or solid");
  }
 
  // auto fe_h = getFEEntityHandle();
  // CHKERR getPtrFE()->mField.get_moab().tag_set_data(th_n, &fe_h, 1,
  //                                                   &t_normal(0));
 
  // calculate  penalty
  const double s = getMeasure() / (areaMap[0] + areaMap[1]);
  const double p = pow(penalty, order) * s;
 
  // get number of integration points on face
  const size_t nb_integration_pts = getGaussPts().size2();
 
  // when boundary edge/face is evaluated, and 2 when is skeleton edge/face.
  const double beta = 1. / (in_the_loop + 1);
 
  auto integrate = [&](auto sense_row, auto &row_ind, auto &row, auto &row_diff,
                       auto sense_col, auto &col_ind, auto &col, auto &col_diff,
                       auto &d_row, auto &d_col, auto &d_vis_row,
                       auto &d_vis_col) {
    MoFEMFunctionBeginHot;
 
    // number of dofs, for homogenous approximation this value is
    // constant.
    const auto nb_rows = row_ind.size();
    const auto nb_cols = col_ind.size();
 
    const auto nb_row_base_functions = row.size2();
 
    auto t_D_row = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(d_row);
    auto t_D_col = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(d_col);
 
    auto t_vis_D_row =
        getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(d_vis_row);
    auto t_vis_D_col =
        getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(d_vis_col);
 
    if (nb_cols) {
 
      // resize local element matrix
      locMatK.resize(nb_rows, nb_cols, false);
      locMatK.clear();
      locMatC.resize(nb_rows, nb_cols, false);
      locMatC.clear();
      locMatM.resize(nb_rows, nb_cols, false);
      locMatM.clear();
 
      // get base functions, and integration weights
      auto t_row_base = get_ntensor(row);
      auto t_diff_row_base = get_diff_ntensor(row_diff);
 
      constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
      FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_P;
      t_P(i, j) = t_normal(i) * t_normal(j);
 
      auto t_w = getFTensor0IntegrationWeight();
 
      // iterate integration points on face/edge
      for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
 
        // t_w is integration weight, and measure is element area, or
        // volume, depending if problem is in 2d/3d.
        const double alpha = getMeasure() * t_w;
 
        // iterate rows
        size_t rr = 0;
        for (; rr != nb_rows / SPACE_DIM; ++rr) {
 
          auto get_vn_plus = [&](auto &t_D) {
            FTensor::Tensor3<double, SPACE_DIM, SPACE_DIM, SPACE_DIM> t_vn_plus;
            t_vn_plus(i, j, k) =
                (beta * phi / p) *
                (t_kd(i, m) * t_D(m, j, k, l) * t_diff_row_base(l));
            return t_vn_plus;
          };
 
          auto get_vn = [&]() {
            FTensor::Tensor3<double, SPACE_DIM, SPACE_DIM, SPACE_DIM> t_vn;
            t_vn(i, j, k) = t_normal(j) * (t_kd(i, k) * t_row_base * sense_row);
            return t_vn;
          };
 
          auto t_vn_plus = get_vn_plus(t_D_row);
          auto t_vn = get_vn();
          auto t_vn_lambda_plus = get_vn_plus(t_vis_D_row);
 
          // get base functions on columns
          auto t_col_base = get_ntensor(col, gg, 0);
          auto t_diff_col_base = get_diff_ntensor(col_diff, gg, 0);
 
          auto t_mat_K = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
              locMatK, SPACE_DIM * rr);
          auto t_mat_C = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
              locMatC, SPACE_DIM * rr);
          auto t_mat_M = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
              locMatM, SPACE_DIM * rr);
 
          // iterate columns
          for (size_t cc = 0; cc != nb_cols / SPACE_DIM; ++cc) {
 
            auto get_un_plus = [&](auto &t_D) {
              FTensor::Tensor3<double, SPACE_DIM, SPACE_DIM, SPACE_DIM>
                  t_un_plus;
              t_un_plus(i, j, k) =
                  beta * ((t_P(i, m) * t_D(m, j, k, l)) * t_diff_col_base(l));
              return t_un_plus;
            };
 
            auto get_un = [&]() {
              FTensor::Tensor3<double, SPACE_DIM, SPACE_DIM, SPACE_DIM> t_un;
              t_un(i, j, k) =
                  (t_normal(j) * (t_P(i, k) * t_col_base * sense_col));
              return t_un;
            };
 
            auto t_un_plus = get_un_plus(t_D_col);
            auto t_un = get_un();
            auto t_un_lambda_plus = get_un_plus(t_vis_D_col);
 
            // K
            t_mat_K(i, j) -=
                alpha * (t_vn(m, n, i) - t_vn_plus(m, n, i)) *
                ((-p) * (t_un(m, n, j) - (t_un_plus(m, n, j) / p)));
            t_mat_K(i, j) -= alpha * (t_vn_plus(m, n, i) * t_un_plus(m, n, j));
 
            // C1
            t_mat_C(i, j) -=
                alpha * (-t_vn_lambda_plus(m, n, i)) *
                ((-p) * (t_un(m, n, j) - (t_un_plus(m, n, j) / p)));
            t_mat_C(i, j) -=
                alpha * (t_vn_lambda_plus(m, n, i) * t_un_plus(m, n, j));
 
            // C2
            t_mat_C(i, j) -= alpha * (t_vn(m, n, i) - t_vn_plus(m, n, i)) *
                             ((-p) * (t_un_lambda_plus(m, n, j) / (-p)));
            t_mat_C(i, j) -=
                alpha * (t_vn_plus(m, n, i) * t_un_lambda_plus(m, n, j));
 
            // M
            t_mat_M(i, j) -= alpha * (-t_vn_lambda_plus(m, n, i)) *
                             ((-p) * (t_un_lambda_plus(m, n, j) / (-p)));
            t_mat_M(i, j) -=
                alpha * (t_vn_lambda_plus(m, n, i) * t_un_lambda_plus(m, n, j));
 
            // move to next column base and element of matrix
            ++t_col_base;
            ++t_diff_col_base;
            ++t_mat_K;
            ++t_mat_C;
            ++t_mat_M;
          }
 
          // move to next row base
          ++t_row_base;
          ++t_diff_row_base;
        }
 
        for (; rr < nb_row_base_functions; ++rr) {
          ++t_row_base;
          ++t_diff_row_base;
        }
 
        ++t_w;
      }
 
      // assemble system
 
      CHKERR ::MatSetValues(K, nb_rows, &*row_ind.begin(), col_ind.size(),
                            &*col_ind.begin(), &*locMatK.data().begin(),
                            ADD_VALUES);
      CHKERR ::MatSetValues(C, nb_rows, &*row_ind.begin(), col_ind.size(),
                            &*col_ind.begin(), &*locMatC.data().begin(),
                            ADD_VALUES);
      CHKERR ::MatSetValues(M, nb_rows, &*row_ind.begin(), col_ind.size(),
                            &*col_ind.begin(), &*locMatM.data().begin(),
                            ADD_VALUES);
    }
 
    MoFEMFunctionReturnHot(0);
  };
 
  // iterate the sides rows
  for (auto s0 : {LEFT_SIDE, RIGHT_SIDE}) {
    const auto sense_row = senseMap[s0];
 
    for (auto x0 = 0; x0 != indices_side_map[s0].size(); ++x0) {
 
      for (auto s1 : {LEFT_SIDE, RIGHT_SIDE}) {
        const auto sense_col = senseMap[s1];
 
        for (auto x1 = 0; x1 != indices_side_map[s1].size(); ++x1) {
 
          CHKERR integrate(sense_row, indices_side_map[s0][x0],
                           base_side_map[s0][x0], diff_base_side_map[s0][x0],
 
                           sense_col, indices_side_map[s1][x1],
                           base_side_map[s1][x1], diff_base_side_map[s1][x1],
 
                           *d_ptr[s0], *d_ptr[s1], *d_vis_ptr[s0],
                           *d_vis_ptr[s1]
 
          );
        }
      }
    }
  }
 
  MoFEMFunctionReturn(0);
}