RotatingFrameOperators.cpp

Go to the documentation of this file.

/** \file RotatingFrameOperators.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 <PlasticOperators.hpp>
#include <RotatingFrameOperators.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 OpRotatingFrameTools {
 
OpRotatingFrameRhs::OpRotatingFrameRhs(
    const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
    : DomainEleOpAssembly(field_name, field_name, DomainEleOp::OPROW),
      commonDataPtr(common_data_ptr) {}
 
MoFEMErrorCode OpRotatingFrameRhs::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_diff_base = data.getFTensor1DiffN<3>();
    auto t_grad = getFTensor2FromMat<3, 3>(*(commonDataPtr->mGradPtr));
    auto t_coords = getFTensor1CoordsAtGaussPts();
 
    for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
 
      Tensor1<double, 3> omg;
      omg(i) = (*cache).Omega(i, j) * t_coords(j);
 
      double alpha = getMeasure() * t_w * (*cache).density;
      auto t_nf = getFTensor1FromArray<3, 3>(locF);
 
      size_t bb = 0;
      for (; bb != nb_dofs / 3; ++bb) {
 
        t_nf(i) -= alpha * (t_grad(i, j) + kronecker_delta(i, j)) * omg(j) *
                   t_diff_base(k) * omg(k);
 
        ++t_diff_base;
        ++t_nf;
      }
      for (; bb < nb_base_functions; ++bb)
        ++t_diff_base;
 
      ++t_coords;
      ++t_grad;
      ++t_w;
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
OpRotatingFrameLhs::OpRotatingFrameLhs(
    const std::string row_field_name, const std::string col_field_name,
    boost::shared_ptr<CommonData> common_data_ptr)
    : DomainEleOpAssembly(row_field_name, col_field_name,
                          DomainEleOp::OPROWCOL),
      commonDataPtr(common_data_ptr) {
  sYmm = false;
}
 
//! [Stiffness]
MoFEMErrorCode OpRotatingFrameLhs::iNtegrate(EntData &row_data,
                                             EntData &col_data) {
  MoFEMFunctionBegin;
 
  const size_t nb_row_dofs = row_data.getIndices().size();
  const size_t nb_col_dofs = col_data.getIndices().size();
 
  if (nb_row_dofs && nb_col_dofs) {
 
    const size_t nb_integration_pts = row_data.getN().size1();
    const size_t nb_row_base_funcs = row_data.getN().size2();
    auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_coords = getFTensor1CoordsAtGaussPts();
    for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
      double alpha = getMeasure() * t_w * (*cache).density;
 
      Tensor1<double, 3> omg;
      omg(i) = (*cache).Omega(i, j) * t_coords(j);
 
      size_t rr = 0;
      for (; rr != nb_row_dofs / 3; ++rr) {
 
        auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
 
        auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
 
        for (size_t cc = 0; cc != nb_col_dofs / 3; ++cc) {
          t_mat(i, k) -= alpha * t_row_diff_base(j) * omg(j) *
                         t_col_diff_base(l) * omg(l) * kronecker_delta(i, k);
 
          ++t_col_diff_base;
          ++t_mat;
        }
 
        ++t_row_diff_base;
      }
      for (; rr != nb_row_base_funcs; ++rr)
        ++t_row_diff_base;
 
      ++t_coords;
      ++t_w;
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
OpRotatingFrameBoundaryRhs::OpRotatingFrameBoundaryRhs(
    const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr,
    boost::shared_ptr<DomainSideEle> side_op_fe)
    : BoundaryEleOpAssembly(field_name, field_name, OPROW),
      commonDataPtr(common_data_ptr), sideOpFe(side_op_fe) {}
 
MoFEMErrorCode OpRotatingFrameBoundaryRhs::iNtegrate(EntData &data) {
  MoFEMFunctionBegin;
 
  const size_t nb_gauss_pts = getGaussPts().size2();
  const size_t nb_dofs = data.getIndices().size();
 
  if (nb_dofs) {
 
    auto t_w = getFTensor0IntegrationWeight();
    auto t_coords = getFTensor1CoordsAtGaussPts();
    CHKERR loopSideVolumes("dFE", *sideOpFe);
    auto t_grad = getFTensor2FromMat<3, 3>(*(commonDataPtr->mGradPtr));
    size_t nb_base_functions = data.getN().size2();
    auto t_base = data.getFTensor0N();
    auto t_normal = getFTensor1Normal();
    const double ls = sqrt(t_normal(i) * t_normal(i));
    t_normal(i) /= ls;
 
    for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
 
      auto t_nf = getFTensor1FromArray<3, 3>(locF);
 
      if (getNormalsAtGaussPts().size1() == nb_gauss_pts) {
        VectorDouble n = getNormalsAtGaussPts(gg);
        auto t_n = getFTensor1FromPtr<3>(&*n.data().begin());
        t_normal(i) = t_n(i) / sqrt(t_n(j) * t_n(j));
      }
 
      const double alpha = t_w * getMeasure() * (*cache).density;
      Tensor1<double, 3> omg;
      omg(i) = (*cache).Omega(i, j) * t_coords(j);
 
      size_t bb = 0;
      for (; bb != nb_dofs / 3; ++bb) {
        // cerr << test << endl;
        t_nf(i) += alpha * t_base *
                   ((t_grad(i, j) + kronecker_delta(i, j)) * omg(j)) *
                   (omg(k) * t_normal(k));
        // alpha * t_base * (t_grad(i, j) * omg(j));
 
        ++t_nf;
        ++t_base;
      }
      for (; bb < nb_base_functions; ++bb)
        ++t_base;
 
      ++t_grad;
      ++t_coords;
      ++t_w;
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
OpRotatingFrameBoundaryLhs::OpRotatingFrameBoundaryLhs(
    const std::string row_field_name, const std::string col_field_name,
    boost::shared_ptr<CommonData> common_data_ptr,
    boost::shared_ptr<DomainSideEle> side_op_fe)
    : BoundaryEleOp(row_field_name, col_field_name, OPROWCOL),
      commonDataPtr(common_data_ptr), sideOpFe(side_op_fe) {
  sYmm = false;
}
 
MoFEMErrorCode OpRotatingFrameBoundaryLhs::doWork(int row_side, int col_side,
                                                  EntityType row_type,
                                                  EntityType col_type,
                                                  EntData &row_data,
                                                  EntData &col_data) {
  MoFEMFunctionBegin;
 
  const size_t nb_gauss_pts = getGaussPts().size2();
  const size_t row_nb_dofs = row_data.getIndices().size();
  size_t col_nb_dofs = col_data.getIndices().size();
 
  if (row_nb_dofs && col_nb_dofs) {
 
    if (col_type == MBVERTEX) {
      commonDataPtr->faceRowData = &row_data;
      commonDataPtr->mMeasure = getMeasure();
      CHKERR loopSideVolumes("dFE", *sideOpFe);
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
OpVolumeSideAssembleRowData::OpVolumeSideAssembleRowData(
    const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
    : DomainSideEleOpAssembly(field_name, field_name, DomainEleOp::OPROWCOL),
      commonDataPtr(common_data_ptr) {
  sYmm = false;
}
 
MoFEMErrorCode OpVolumeSideAssembleRowData::iNtegrate(EntData &row_data,
                                                      EntData &col_data) {
  return 0;
}
 
MoFEMErrorCode OpVolumeSideAssembleRowData::doWork(
    int row_side, int col_side, EntityType row_type, EntityType col_type,
    EntitiesFieldData::EntData &row_data,
    EntitiesFieldData::EntData &col_data) {
  MoFEMFunctionBegin;
  if (row_type != MBVERTEX)
    MoFEMFunctionReturnHot(0);
 
  const size_t nb_gauss_pts = getGaussPts().size2();
  const size_t row_nb_dofs = commonDataPtr->faceRowData->getIndices().size();
  const size_t col_nb_dofs = col_data.getIndices().size();
 
  if (row_nb_dofs && col_nb_dofs) {
 
    auto t_coords = getFTensor1CoordsAtGaussPts();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_row_base = commonDataPtr->faceRowData->getFTensor0N();
 
    auto t_row_diff_base = commonDataPtr->faceRowData->getFTensor1DiffN<3>();
    size_t nb_face_functions = commonDataPtr->faceRowData->getN().size2();
    auto t_normal = getFTensor1Normal();
    const double ls = sqrt(t_normal(i) * t_normal(i));
    t_normal(i) /= ls;
 
    locMat.resize(row_nb_dofs, col_nb_dofs, false);
    locMat.clear();
 
    for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
      // we are actually integrating on a triangle here, therefore 0.5
      // const double alpha_test = t_w * getMeasure() * (*cache).density;
      const double alpha = t_w * commonDataPtr->mMeasure * (*cache).density;
      // const double alpha = t_w * 0.5 * (*cache).density;
      Tensor1<double, 3> omg;
      omg(i) = (*cache).Omega(i, j) * t_coords(j);
 
      if (getNormalsAtGaussPts().size1() == nb_gauss_pts) {
        VectorDouble n = getNormalsAtGaussPts(gg);
        auto t_n = getFTensor1FromPtr<3>(&*n.data().begin());
        t_normal(i) = t_n(i) / sqrt(t_n(j) * t_n(j));
      }
 
      size_t rr = 0;
      for (; rr != row_nb_dofs / 3; ++rr) {
        auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
 
        const double beta = alpha * t_row_base;
        auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
 
        for (size_t cc = 0; cc != col_nb_dofs / 3; ++cc) {
 
          t_mat(i, k) += beta * kronecker_delta(i, k) *
                         (omg(j) * t_col_diff_base(j)) * (omg(l) * t_normal(l));
 
          ++t_col_diff_base;
          ++t_mat;
        }
 
        ++t_row_base;
      }
      for (; rr < nb_face_functions; ++rr)
        ++t_row_base;
 
      ++t_coords;
      ++t_w;
    }
    // CHKERR iNtegrate(row_data, col_data);
    CHKERR aSsemble(*commonDataPtr->faceRowData, col_data, false);
  }
  MoFEMFunctionReturn(0);
}
 
OpCalculatePlasticConvRotatingFrame::OpCalculatePlasticConvRotatingFrame(
    const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
    : DomainEleOp(field_name, field_name, DomainEleOp::OPROW),
      commonDataPtr(common_data_ptr) {}
 
MoFEMErrorCode OpCalculatePlasticConvRotatingFrame::doWork(int side,
                                                           EntityType type,
                                                           EntData &data) {
  MoFEMFunctionBegin;
  const size_t nb_dofs = data.getIndices().size();
  if (nb_dofs) {
    const size_t nb_integration_pts = data.getN().size1();
 
    auto t_plastic_strain_dot =
        getFTensor2SymmetricFromMat<3>(*(commonDataPtr->plasticStrainDotPtr));
    auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
    auto t_coords = getFTensor1CoordsAtGaussPts();
    auto t_grad_tau =
        getFTensor1FromMat<3>(*(commonDataPtr->plasticGradTauPtr));
 
    // OpCalculateTensor2SymmetricFieldGradient
    auto t_grad_plastic_strain =
        getFTensor3DgFromMat<3, 3>(*(commonDataPtr->plasticGradStrainPtr));
    commonDataPtr->guidingVelocityPtr->resize(3, nb_integration_pts, false);
    commonDataPtr->guidingVelocityPtr->clear();
    auto t_omega = getFTensor1FromMat<3>(*commonDataPtr->guidingVelocityPtr);
 
    for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
 
      t_omega(i) = (*cache).Omega(i, j) * t_coords(j);
      t_tau_dot += t_grad_tau(i) * t_omega(i);
      t_plastic_strain_dot(i, j) += t_grad_plastic_strain(i, j, k) * t_omega(k);
 
      ++t_coords;
      ++t_omega;
      ++t_grad_plastic_strain;
      ++t_plastic_strain_dot;
      ++t_tau_dot;
      ++t_grad_tau;
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
OpCalculateJumpOnSkeleton::OpCalculateJumpOnSkeleton(
    const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr,
    boost::shared_ptr<DomainSideEle> side_op_fe)
    : SkeletonEleOp(field_name, field_name, SkeletonEleOp::OPROW),
      commonDataPtr(common_data_ptr), sideOpFe(side_op_fe) {
  std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
  doEntities[MBTRI] = doEntities[MBQUAD] = true;
}
 
MoFEMErrorCode OpCalculateJumpOnSkeleton::doWork(int side, EntityType type,
                                                 EntData &data) {
  MoFEMFunctionBegin;
  // if (skipEnts.find(getFEEntityHandle()) != skipEnts.end())
  //   MoFEMFunctionReturnHot(0);
  if ((type == MBTRI || type == MBQUAD) && side == 0) {
    const size_t nb_dofs = data.getIndices().size();
    // if (nb_dofs == 0)
    //   MoFEMFunctionReturnHot(0);
 
    const size_t nb_integration_pts = data.getN().size1();
    const size_t nb_base_functions = data.getN().size2();
    commonDataPtr->plasticStrainSideMap.clear();
    CHKERR loopSideVolumes("dFE", *sideOpFe);
    const int check_size = commonDataPtr->plasticStrainSideMap.size();
 
    if (commonDataPtr->plasticStrainSideMap.size() == 1)
      if (static_cast<int>(
              commonDataPtr->plasticStrainSideMap.begin()->second.size2()) !=
          nb_integration_pts)
        SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                "wrong number of integration points");
 
    auto &plastic_mat_l = commonDataPtr->plasticStrainSideMap.at(-1);
    auto &plastic_mat_r = commonDataPtr->plasticStrainSideMap.at(1);
    auto t_plastic_strain_l = getFTensor2SymmetricFromMat<3>(plastic_mat_l);
    auto t_plastic_strain_r = getFTensor2SymmetricFromMat<3>(plastic_mat_r);
    auto &tau_vec_l = commonDataPtr->plasticTauSideMap.at(-1);
    auto &tau_vec_r = commonDataPtr->plasticTauSideMap.at(1);
    auto t_tau_l = getFTensor0FromVec(tau_vec_l);
    auto t_tau_r = getFTensor0FromVec(tau_vec_r);
 
    auto t_w = getFTensor0IntegrationWeight();
    auto t_base = data.getFTensor0N();
    auto t_coords = getFTensor1CoordsAtGaussPts();
    commonDataPtr->plasticTauJumpPtr->resize(nb_integration_pts, false);
    commonDataPtr->plasticStrainJumpPtr->resize(6, nb_integration_pts, false);
    auto t_plastic_strain_jump =
        getFTensor2SymmetricFromMat<3>(*commonDataPtr->plasticStrainJumpPtr);
    auto t_tau_jump = getFTensor0FromVec(*commonDataPtr->plasticTauJumpPtr);
 
    for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
      double alpha = getMeasure() * t_w;
 
      // Tensor2_symmetric<PackPtr<double *, 6>, 3> t_nf{&nf[0], &nf[1],
      // &nf[2],
      //                                                 &nf[3], &nf[4],
      //                                                 &nf[5]};
      // this is just work in progress
      t_plastic_strain_jump(i, j) =
          t_plastic_strain_l(i, j) - t_plastic_strain_r(i, j);
      t_tau_jump = t_tau_l - t_tau_r;
 
      // size_t bb = 0;
      // for (; bb != nb_dofs / 6; ++bb) {
      //   t_nf(i, j) += alpha * t_base;
      //   ++t_base;
      //   ++t_nf;
      // }
 
      // for (; bb < nb_base_functions; ++bb)
      //   ++t_base;
 
      ++t_plastic_strain_jump;
      ++t_tau_jump;
      ++t_plastic_strain_l;
      ++t_plastic_strain_r;
      ++t_tau_l;
      ++t_tau_r;
      ++t_w;
    }
 
    // CHKERR VecSetValues(getTSf(), data, nf.data(), ADD_VALUES);
  }
 
  MoFEMFunctionReturn(0);
}
 
OpVolumeSideCalculateEP::OpVolumeSideCalculateEP(
    const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
    : DomainSideEleOp(field_name, field_name, DomainEleOp::OPROW),
      commonDataPtr(common_data_ptr) {
  std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
  doEntities[MBTET] = true;
  doEntities[MBHEX] = true;
}
 
MoFEMErrorCode OpVolumeSideCalculateEP::doWork(int side, EntityType type,
                                               EntData &data) {
  MoFEMFunctionBegin;
  const int nb_gauss_pts = getGaussPts().size2();
  const int nb_base_functions = data.getN().size2();
  // commonDataPtr->mStressPtr->resize(6, nb_gauss_pts);
  // auto &t_D = commonDataPtr->tD;
  // auto t_strain =
  // getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStrainPtr)); auto t_stress
  // = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStressPtr));
  // FIXME: this is wrong, resize only on zero type
  const size_t nb_dofs = data.getIndices().size();
  // data.getFieldData().size() == 0
  // if (side == getFaceSideNumber()) {
  if (nb_dofs) {
 
    // for testing //FIXME:
    MatrixDouble diff = getCoordsAtGaussPts() - getFaceCoordsAtGaussPts();
    const double eps = 1e-12;
    if (norm_inf(diff) > eps)
      SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
              "coordinates at integration pts are different");
 
    const size_t nb_gauss_pts = data.getN().size1();
 
    auto base_function = data.getFTensor0N();
    // MatrixDouble *ep_elem_data = nullptr;
    int test_nb_in_loop = getFEMethod()->nInTheLoop;
    int test_face_sense = getFaceSense();
    MatrixDouble &mat_ep = commonDataPtr->plasticStrainSideMap[getFaceSense()];
    mat_ep.resize(6, nb_gauss_pts, false);
    mat_ep.clear();
 
    // if (getFEMethod()->nInTheLoop == 0)
    //   ep_elem_data = &(commonDataPtr->plasticStrainLeft);
    // else
    //   ep_elem_data = &(commonDataPtr->plasticStrainRight);
 
    // mat_ep.resize(6, nb_gauss_pts, false);
    // mat_ep.clear();
    auto values_at_gauss_pts = getFTensor2SymmetricFromMat<3>(mat_ep);
    for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
      auto field_data = data.getFTensor2SymmetricFieldData<3>();
      size_t bb = 0;
      for (; bb != nb_dofs / 6; ++bb) {
        values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
        ++field_data;
        ++base_function;
      }
      for (; bb != nb_base_functions; ++bb)
        ++base_function;
      ++values_at_gauss_pts;
    }
  }
  MoFEMFunctionReturn(0);
}
 
OpVolumeSideCalculateTAU::OpVolumeSideCalculateTAU(
    const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
    : DomainSideEleOp(field_name, field_name, DomainEleOp::OPROW),
      commonDataPtr(common_data_ptr) {
  std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
  doEntities[MBTET] = true;
  doEntities[MBHEX] = true;
}
 
MoFEMErrorCode OpVolumeSideCalculateTAU::doWork(int side, EntityType type,
                                                EntData &data) {
  MoFEMFunctionBegin;
 
  // FIXME: this is wrong, resize only on zero type
  const size_t nb_dofs = data.getIndices().size();
  // data.getFieldData().size()
  if (nb_dofs) {
    int nb_gauss_pts = data.getN().size1();
    VectorDouble &vec_tau = commonDataPtr->plasticTauSideMap[getFaceSense()];
    vec_tau.resize(nb_gauss_pts, false);
    vec_tau.clear();
 
    for (int gg = 0; gg != nb_gauss_pts; gg++) {
      vec_tau(gg) = inner_prod(trans(data.getN(gg)), data.getFieldData());
    }
  }
  MoFEMFunctionReturn(0);
}
 
OpPostProcPlasticSkeleton::OpPostProcPlasticSkeleton(
    const std::string field_name, moab::Interface &post_proc_mesh,
    std::vector<EntityHandle> &map_gauss_pts,
    boost::shared_ptr<CommonData> common_data_ptr)
    : SkeletonEleOp(field_name, field_name, SkeletonEleOp::OPROW),
      postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
      commonDataPtr(common_data_ptr) {
  // Operator is only executed for vertices
  std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
  doEntities[MBTRI] = doEntities[MBQUAD] = true;
}
 
//! [Postprocessing]
MoFEMErrorCode OpPostProcPlasticSkeleton::doWork(int side, EntityType type,
                                                 EntData &data) {
  MoFEMFunctionBegin;
  // if (skinEnts.find(getFEEntityHandle()) != skinEnts.end())
  //   MoFEMFunctionReturnHot(0);
  const size_t nb_gauss_pts = getGaussPts().size2();
  std::array<double, 9> def;
  std::fill(def.begin(), def.end(), 0);
 
  auto get_tag = [&](const std::string name, size_t size) {
    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_3d = [&](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_vector = [&](auto &t) -> MatrixDouble3by3 & {
    mat.clear();
    for (size_t r = 0; r != 3; ++r)
      mat(0, r) = t(r);
    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());
  };
 
  auto th_rot_vel = get_tag("ROT_VELOCITY", 3);
  auto th_tau_jump = get_tag("TAU_JUMP", 1);
  auto th_plastic_jump = get_tag("PLASTIC_STRAIN_JUMP", 9);
  auto t_coords = getFTensor1CoordsAtGaussPts();
  auto t_ep_jump =
      getFTensor2SymmetricFromMat<3>(*(commonDataPtr->plasticStrainJumpPtr));
  auto t_tau_jump = getFTensor0FromVec(*(commonDataPtr->plasticTauJumpPtr));
 
  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
 
    Tensor1<double, 3> velocity;
    velocity(i) = (*cache).Omega(i, j) * t_coords(j);
 
    CHKERR set_tag(th_tau_jump, gg, set_scalar(t_tau_jump));
    CHKERR set_tag(th_plastic_jump, gg, set_matrix_3d(t_ep_jump));
    CHKERR set_tag(th_rot_vel, gg, set_vector(velocity));
    ++t_coords;
    ++t_ep_jump;
    ++t_tau_jump;
  }
 
  MoFEMFunctionReturn(0);
}
 
} // namespace OpRotatingFrameTools