CrackFrontElement.cpp

Go to the documentation of this file.

/** \file CrackFrontElement.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 <moab/SpatialLocator.hpp>
#include <moab/ElemEvaluator.hpp>
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <cholesky.hpp>
 
#include <BasicFiniteElements.hpp>
#include <Mortar.hpp>
#include <Hooke.hpp>
 
 
#include <AnalyticalFun.hpp>
 
 
// TODO implementation in ComplexConstArea is obsolete, not follows naming
// convention and code is unnecessarily complicated and in places difficult to
// follow.
#include <ComplexConstArea.hpp>
#include <ConstantArea.hpp>
 
extern "C" {
#include <phg-quadrule/quad.h>
}
 
#include <NeoHookean.hpp>
#include <Hooke.hpp>
 
#include <CrackFrontElement.hpp>
 
namespace po = boost::program_options;
 
namespace FractureMechanics {
 
static std::map<long int, MatrixDouble> mapRefCoords;
 
template <>
template <>
int CrackFrontSingularBase<NonlinearElasticElement::MyVolumeFE,
                           VolumeElementForcesAndSourcesCore>::
    getRuleImpl(TwoType<NonlinearElasticElement::MyVolumeFE,
                        VolumeElementForcesAndSourcesCore>,
                int order) {
  return -1;
}
 
template <>
template <>
MoFEMErrorCode CrackFrontSingularBase<NonlinearElasticElement::MyVolumeFE,
                                      VolumeElementForcesAndSourcesCore>::
    setGaussPtsImpl(TwoType<NonlinearElasticElement::MyVolumeFE,
                            VolumeElementForcesAndSourcesCore>,
                    int order) {
 
  MoFEMFunctionBeginHot;
 
  auto set_gauss_pts = [this]() {
    MoFEMFunctionBegin;
    gaussPts.swap(refGaussCoords);
    const size_t nb_gauss_pts = gaussPts.size2();
    dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4);
    double *shape_ptr =
        &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
    CHKERR ShapeMBTET(shape_ptr, &gaussPts(0, 0), &gaussPts(1, 0),
                      &gaussPts(2, 0), nb_gauss_pts);
    MoFEMFunctionReturn(0);
  };
 
  std::bitset<4> singular_nodes_ids;
  if (refineIntegration == PETSC_TRUE && singularElement) {
 
    for (int nn = 0; nn != 4; nn++) {
      if (singularNodes[nn])
        singular_nodes_ids.set(nn);
    }
 
    auto it_map_ref_coords = mapRefCoords.find(singular_nodes_ids.to_ulong());
    if (it_map_ref_coords != mapRefCoords.end()) {
 
      refGaussCoords = it_map_ref_coords->second;
      CHKERR set_gauss_pts();
 
      MoFEMFunctionReturnHot(0);
    }
  }
 
  order = order ? order : 1;
  int rule = (singularElement) ? 2 * (std::max(order, 3) - 1) : 2 * (order - 1);
 
  if (rule < QUAD_3D_TABLE_SIZE) {
    if (QUAD_3D_TABLE[rule]->dim != 3) {
      SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "wrong dimension");
    }
    if (QUAD_3D_TABLE[rule]->order < rule) {
      SETERRQ2(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
               "wrong order %d != %d", QUAD_3D_TABLE[rule]->order, rule);
    }
    const size_t nb_gauss_pts = QUAD_3D_TABLE[rule]->npoints;
    gaussPts.resize(4, nb_gauss_pts, false);
    cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[1], 4, &gaussPts(0, 0),
                1);
    cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[2], 4, &gaussPts(1, 0),
                1);
    cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[3], 4, &gaussPts(2, 0),
                1);
    cblas_dcopy(nb_gauss_pts, QUAD_3D_TABLE[rule]->weights, 1, &gaussPts(3, 0),
                1);
    dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4,
                                                           false);
    double *shape_ptr =
        &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
    cblas_dcopy(4 * nb_gauss_pts, QUAD_3D_TABLE[rule]->points, 1, shape_ptr, 1);
  } else {
    SETERRQ2(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
             "rule > quadrature order %d < %d", rule, QUAD_3D_TABLE_SIZE);
  }
 
  // Refine quadrature by splitting edges adjacent to crack front. On sub
  // triangle place quadrature points.
  if (!singularElement)
    MoFEMFunctionReturnHot(0);
 
  // refine mesh on refine elements
  if (refineIntegration == PETSC_TRUE) {
 
    const int max_level = refinementLevels;
    EntityHandle tet;
 
    moab::Core moab_ref;
    double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
    EntityHandle nodes[4];
    for (int nn = 0; nn != 4; nn++)
      CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
    CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
    MoFEM::CoreTmp<-1> mofem_ref_core(moab_ref, PETSC_COMM_SELF, -2);
    MoFEM::Interface &m_field_ref = mofem_ref_core;
    CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevelByDim(
        0, 3, BitRefLevel().set(0));
 
    Range singular_nodes;
    for (int nn = 0; nn != 4; nn++) {
      if (singularNodes[nn]) {
        EntityHandle ent;
        CHKERR moab_ref.side_element(tet, 0, nn, ent);
        singular_nodes.insert(ent);
      }
    }
 
    EntityHandle meshset;
    if (singular_nodes.size() > 1) {
      CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
      for (int ee = 0; ee != 6; ee++) {
        if (singularEdges[ee]) {
          EntityHandle ent;
          CHKERR moab_ref.side_element(tet, 1, ee, ent);
          CHKERR moab_ref.add_entities(meshset, &ent, 1);
        }
      }
    }
 
    MeshRefinement *m_ref;
    CHKERR m_field_ref.getInterface(m_ref);
    for (int ll = 0; ll != max_level; ll++) {
      // PetscPrintf(T::mField.get_comm(),"Refine Level %d\n",ll);
      Range edges;
      CHKERR m_field_ref.getInterface<BitRefManager>()
          ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
                                         BitRefLevel().set(), MBEDGE, edges);
      Range ref_edges;
      CHKERR moab_ref.get_adjacencies(singular_nodes, 1, true, ref_edges,
                                      moab::Interface::UNION);
      ref_edges = intersect(ref_edges, edges);
      if (singular_nodes.size() > 1) {
        Range ents;
        CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents, true);
        ref_edges = intersect(ref_edges, ents);
      }
      Range tets;
      CHKERR m_field_ref.getInterface<BitRefManager>()
          ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
                                         BitRefLevel().set(), MBTET, tets);
      // refine mesh
      CHKERR m_ref->addVerticesInTheMiddleOfEdges(
          ref_edges, BitRefLevel().set(ll + 1));
      CHKERR m_ref->refineTets(tets, BitRefLevel().set(ll + 1));
      if (singular_nodes.size() > 1) {
        CHKERR m_field_ref.getInterface<BitRefManager>()
            ->updateMeshsetByEntitiesChildren(
                meshset, BitRefLevel().set(ll + 1), meshset, MBEDGE, true);
      }
    }
 
    Range tets;
    CHKERR m_field_ref.getInterface<BitRefManager>()
        ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(max_level),
                                       BitRefLevel().set(), MBTET, tets);
    refCoords.resize(tets.size(), 12, false);
    int tt = 0;
    for (Range::iterator tit = tets.begin(); tit != tets.end(); tit++, tt++) {
      int num_nodes;
      const EntityHandle *conn;
      CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
      CHKERR moab_ref.get_coords(conn, num_nodes, &refCoords(tt, 0));
    }
 
    if (0) {
      EntityHandle out_meshset;
      CHKERR moab_ref.create_meshset(MESHSET_SET, out_meshset);
      CHKERR moab_ref.add_entities(out_meshset, tets);
      std::string file_name =
          "ref_tet_" + boost::lexical_cast<std::string>(nInTheLoop) + ".vtk";
      CHKERR moab_ref.write_file(file_name.c_str(), "VTK", "", &out_meshset, 1);
      CHKERR moab_ref.delete_entities(&out_meshset, 1);
    }
 
    // Set integration points
    double diff_n[12];
    CHKERR ShapeDiffMBTET(diff_n);
    const size_t nb_gauss_pts = gaussPts.size2();
    refGaussCoords.resize(4, nb_gauss_pts * refCoords.size1());
    MatrixDouble &shape_n = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
    int gg = 0;
    for (size_t tt = 0; tt != refCoords.size1(); tt++) {
      double *tet_coords = &refCoords(tt, 0);
      double det = ShapeVolumeMBTET(diff_n, tet_coords);
      det *= 6;
      for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
        for (int dd = 0; dd != 3; dd++) {
          refGaussCoords(dd, gg) = shape_n(ggg, 0) * tet_coords[3 * 0 + dd] +
                                   shape_n(ggg, 1) * tet_coords[3 * 1 + dd] +
                                   shape_n(ggg, 2) * tet_coords[3 * 2 + dd] +
                                   shape_n(ggg, 3) * tet_coords[3 * 3 + dd];
        }
        refGaussCoords(3, gg) = gaussPts(3, ggg) * det;
      }
    }
 
    mapRefCoords[singular_nodes_ids.to_ulong()] = refGaussCoords;
    CHKERR set_gauss_pts();
  }
 
  MoFEMFunctionReturnHot(0);
}
 
MoFEMErrorCode OpGetCrackFrontDataGradientAtGaussPts::doWork(
    int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
  MoFEMFunctionBegin;
 
  if (!singularElement) {
    return OpCalculateVectorFieldGradient<3, 3>::doWork(side, type, data);
  }
 
  const int nb_dofs = data.getFieldData().size();
  if (!nb_dofs && type == this->zeroType) {
    this->dataPtr->resize(9, 0, false);
    MoFEMFunctionReturnHot(0);
  }
  if (!nb_dofs) {
    MoFEMFunctionReturnHot(0);
  }
  const int nb_gauss_pts = data.getN().size1();
  const int nb_base_functions = data.getN().size2();
  ublas::matrix<double, ublas::row_major, DoubleAllocator> &mat =
      *this->dataPtr;
  if (type == this->zeroType) {
    mat.resize(9, nb_gauss_pts, false);
    mat.clear();
  }
  auto diff_base_function = data.getFTensor1DiffN<3>();
  auto gradients_at_gauss_pts = getFTensor2FromMat<3, 3>(mat);
  FTensor::Index<'I', 3> I;
  FTensor::Index<'J', 3> J;
  FTensor::Index<'K', 3> K;
  int size = nb_dofs / 3;
  if (nb_dofs % 3) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
  }
  double *inv_jac_ptr = &*invSJac.data().begin();
  FTensor::Tensor2<double *, 3, 3> t_inv_singular_jacobian(
      inv_jac_ptr, &inv_jac_ptr[1], &inv_jac_ptr[2], &inv_jac_ptr[3],
      &inv_jac_ptr[4], &inv_jac_ptr[5], &inv_jac_ptr[6], &inv_jac_ptr[7],
      &inv_jac_ptr[8], 9);
 
  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
    auto field_data = data.getFTensor1FieldData<3>();
    int bb = 0;
    for (; bb < size; ++bb) {
      gradients_at_gauss_pts(I, J) +=
          field_data(I) *
          (t_inv_singular_jacobian(K, J) * diff_base_function(K));
      ++field_data;
      ++diff_base_function;
    }
    // Number of dofs can be smaller than number of Tensor_Dim0 x base functions
    for (; bb != nb_base_functions; ++bb)
      ++diff_base_function;
    ++gradients_at_gauss_pts;
    ++t_inv_singular_jacobian;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
OpGetCrackFrontDataAtGaussPts::doWork(int side, EntityType type,
                                      DataForcesAndSourcesCore::EntData &data) {
  MoFEMFunctionBegin;
 
  if (!singularElement) {
    return NonlinearElasticElement::OpGetDataAtGaussPts::doWork(side, type,
                                                                data);
  }
 
  const int nb_dofs = data.getFieldData().size();
  const int nb_base_functions = data.getN().size2();
  if (nb_dofs == 0) {
    MoFEMFunctionReturnHot(0);
  }
  const int nb_gauss_pts = data.getN().size1();
  const int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
 
  // initialize
  if (type == zeroAtType) {
    valuesAtGaussPts.resize(nb_gauss_pts);
    gradientAtGaussPts.resize(nb_gauss_pts);
    for (int gg = 0; gg != nb_gauss_pts; gg++) {
      valuesAtGaussPts[gg].resize(rank, false);
      gradientAtGaussPts[gg].resize(rank, 3, false);
    }
    for (int gg = 0; gg != nb_gauss_pts; gg++) {
      valuesAtGaussPts[gg].clear();
      gradientAtGaussPts[gg].clear();
    }
  }
 
  auto base_function = data.getFTensor0N();
  auto diff_base_functions = data.getFTensor1DiffN<3>();
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
 
  double *inv_jac_ptr = &*invSJac.data().begin();
  FTensor::Tensor2<double *, 3, 3> t_inv_singular_jacobian(
      inv_jac_ptr, &inv_jac_ptr[1], &inv_jac_ptr[2], &inv_jac_ptr[3],
      &inv_jac_ptr[4], &inv_jac_ptr[5], &inv_jac_ptr[6], &inv_jac_ptr[7],
      &inv_jac_ptr[8], 9);
 
  if (rank == 3) {
 
    for (int gg = 0; gg != nb_gauss_pts; gg++) {
      auto field_data = data.getFTensor1FieldData<3>();
      FTensor::Tensor1<double *, 3> values(&valuesAtGaussPts[gg][0],
                                           &valuesAtGaussPts[gg][1],
                                           &valuesAtGaussPts[gg][2]);
      FTensor::Tensor2<double *, 3, 3> gradient(
          &gradientAtGaussPts[gg](0, 0), &gradientAtGaussPts[gg](0, 1),
          &gradientAtGaussPts[gg](0, 2), &gradientAtGaussPts[gg](1, 0),
          &gradientAtGaussPts[gg](1, 1), &gradientAtGaussPts[gg](1, 2),
          &gradientAtGaussPts[gg](2, 0), &gradientAtGaussPts[gg](2, 1),
          &gradientAtGaussPts[gg](2, 2));
      int bb = 0;
      for (; bb != nb_dofs / 3; bb++) {
        values(i) += base_function * field_data(i);
        gradient(i, j) += field_data(i) * (t_inv_singular_jacobian(k, j) *
                                           diff_base_functions(k));
        ++diff_base_functions;
        ++base_function;
        ++field_data;
      }
      for (; bb != nb_base_functions; bb++) {
        ++diff_base_functions;
        ++base_function;
      }
      ++t_inv_singular_jacobian;
    }
 
  } else {
 
    SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
             "Not implemented for rank  = %d", rank);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
OpPostProcDisplacements::doWork(int side, EntityType type,
                                DataForcesAndSourcesCore::EntData &data) {
  MoFEMFunctionBegin;
  if (!singularElement || type != MBVERTEX) {
    MoFEMFunctionReturnHot(0);
  }
  double def_VAL[3] = {0, 0, 0};
  Tag th_singular_disp;
  CHKERR postProcMesh.tag_get_handle("SINGULAR_DISP", 3, MB_TYPE_DOUBLE,
                                     th_singular_disp,
                                     MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  auto t_H = getFTensor2FromMat<3, 3>(*H);
  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_singular_disp(
      &singularDisp(0, 0), &singularDisp(0, 1), &singularDisp(0, 2));
  FTensor::Tensor1<double, 3> t_singular_current_disp;
  for (unsigned int gg = 0; gg != singularDisp.size1(); gg++) {
    t_singular_current_disp(i) = t_H(i, j) * t_singular_disp(j);
    CHKERR postProcMesh.tag_set_data(th_singular_disp, &mapGaussPts[gg], 1,
                                     &t_singular_current_disp(0));
    ++t_H;
    ++t_singular_disp;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpTransfromSingularBaseFunctions::doWork(
    int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  MoFEMFunctionBegin;
 
  if (!singularElement)
    MoFEMFunctionReturnHot(0);
 
  if (invSJac.size2() != 9) {
    SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
             "It looks that ho inverse of Jacobian is not calculated %d != 9",
             invSJac.size2());
  }
  double *t_inv_jac_ptr = &*invSJac.data().begin();
  FTensor::Tensor2<double *, 3, 3> t_inv_jac(
      t_inv_jac_ptr, &t_inv_jac_ptr[1], &t_inv_jac_ptr[2], &t_inv_jac_ptr[3],
      &t_inv_jac_ptr[4], &t_inv_jac_ptr[5], &t_inv_jac_ptr[6],
      &t_inv_jac_ptr[7], &t_inv_jac_ptr[8], 9);
 
  switch (type) {
  case MBVERTEX:
    if (applyDet) {
      for (unsigned int gg = 0; gg != getGaussPts().size2(); ++gg) {
        getGaussPts()(3, gg) *= detS[gg];
      }
    }
  case MBEDGE:
  case MBTRI:
  case MBTET: {
    unsigned int nb_gauss_pts = data.getN(bAse).size1();
    if (nb_gauss_pts == 0)
      MoFEMFunctionReturnHot(0);
    if (invSJac.size1() != nb_gauss_pts) {
      SETERRQ2(
          PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
          "It looks that ho inverse of Jacobian is not calculated %d != %d",
          invSJac.size1(), nb_gauss_pts);
    }
    unsigned int nb_base_functions = data.getN(bAse).size2();
    if (nb_base_functions == 0)
      MoFEMFunctionReturnHot(0);
    diffNinvJac.resize(nb_gauss_pts, 3 * nb_base_functions, false);
    double *t_inv_n_ptr = &*diffNinvJac.data().begin();
    FTensor::Tensor1<double *, 3> t_inv_diff_n(t_inv_n_ptr, &t_inv_n_ptr[1],
                                               &t_inv_n_ptr[2], 3);
    // Apply transformation to each bAse function ob given element
    FTensor::Tensor1<double *, 3> t_diff_n = FTensor::Tensor1<double *, 3>(
        &data.getDiffN(bAse)(0, 0), &data.getDiffN(bAse)(0, 1),
        &data.getDiffN(bAse)(0, 2), 3);
    for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
      for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
        t_inv_diff_n(i) = t_diff_n(j) * t_inv_jac(j, i);
        ++t_diff_n;
        ++t_inv_diff_n;
      }
      ++t_inv_jac;
    }
    if (type == MBVERTEX) {
      data.getDiffN(bAse).resize(diffNinvJac.size1(), diffNinvJac.size2(),
                                 false);
    }
    data.getDiffN(bAse).data().swap(diffNinvJac.data());
  } break;
  default:
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
  }
 
  MoFEMFunctionReturn(0);
}
 
OpRhsBoneExplicitDerivariveWithHooke::OpRhsBoneExplicitDerivariveWithHooke(
    HookeElement::DataAtIntegrationPts &common_data,
    boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
    boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts, Range tets_in_block,
    const double rho_0, const double bone_n, Range *forces_only_on_entities_row)
    : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
          "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
      commonData(common_data), rhoAtGaussPts(rho_at_gauss_pts),
      rhoGradAtGaussPts(rho_grad_at_gauss_pts), tetsInBlock(tets_in_block),
      rHo0(rho_0), boneN(bone_n) {
  if (forces_only_on_entities_row)
    forcesOnlyOnEntitiesRow = *forces_only_on_entities_row;
}
 
MoFEMErrorCode OpRhsBoneExplicitDerivariveWithHooke::doWork(
    int row_side, EntityType row_type,
    DataForcesAndSourcesCore::EntData &row_data) {
  FTensor::Index<'i', 3> i;
  MoFEMFunctionBeginHot;
  const int nb_dofs = row_data.getIndices().size();
  if (!nb_dofs)
    MoFEMFunctionReturnHot(0);
  if (tetsInBlock.find(getFEEntityHandle()) == tetsInBlock.end()) {
    MoFEMFunctionReturnHot(0);
  }
  // If this is not part of blockset
  const int nb_gauss_pts = row_data.getN().size1();
  nF.resize(nb_dofs, false);
  nF.clear();
  auto t_rho_gradient = getFTensor1FromMat<3>(*rhoGradAtGaussPts);
  auto rho = getFTensor0FromVec(*rhoAtGaussPts);
  auto energy = getFTensor0FromVec(*commonData.energyVec);
  auto det_H = getFTensor0FromVec(*commonData.detHVec);
 
  for (int gg = 0; gg != nb_gauss_pts; gg++) {
    const double w = getVolume() * getGaussPts()(3, gg) * det_H;
 
    const double denergy_drho =
        energy * (boneN / rho); // derivative of energy w.r.t density
    auto t_base = row_data.getFTensor0N(gg, 0);
    FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_nf(&nF[0], &nF[1],
                                                            &nF[2]);
    for (int bb = 0; bb != nb_dofs / 3; bb++) {
      t_nf(i) += t_base * (w * denergy_drho * t_rho_gradient(i));
      ++t_nf;
      ++t_base;
    }
    ++t_rho_gradient;
    ++rho;
    ++energy;
    ++det_H;
  }
  int *indices_ptr = &row_data.getIndices()[0];
  if (!forcesOnlyOnEntitiesRow.empty()) {
    iNdices.resize(nb_dofs, false);
    noalias(iNdices) = row_data.getIndices();
    indices_ptr = &iNdices[0];
    auto &dofs = row_data.getFieldDofs();
    auto dit = dofs.begin();
    for (int ii = 0; dit != dofs.end(); dit++, ii++) {
      if (auto dof = (*dit)) {
        if (forcesOnlyOnEntitiesRow.find(dof->getEnt()) ==
            forcesOnlyOnEntitiesRow.end()) {
          iNdices[ii] = -1;
        }
      }
    }
  }
 
  CHKERR VecSetOption(getFEMethod()->snes_f, VEC_IGNORE_NEGATIVE_INDICES,
                      PETSC_TRUE);
  CHKERR VecSetValues(getFEMethod()->snes_f, nb_dofs, indices_ptr, &nF[0],
                      ADD_VALUES);
  MoFEMFunctionReturnHot(0);
}
 
OpLhsBoneExplicitDerivariveWithHooke_dX::
    OpLhsBoneExplicitDerivariveWithHooke_dX(
        HookeElement::DataAtIntegrationPts &common_data,
        boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
        boost::shared_ptr<MatrixDouble> diff_rho_at_gauss_pts,
        boost::shared_ptr<MatrixDouble> diff_diff_rho_at_gauss_pts,
        MatrixDouble &singular_displ, Range tets_in_block, const double rho_0,
        const double bone_n, Range *forces_only_on_entities_row)
    : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
          "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
          UserDataOperator::OPROWCOL, false), // FIXME:
      commonData(common_data), rhoAtGaussPts(rho_at_gauss_pts),
      diffRhoAtGaussPts(diff_rho_at_gauss_pts),
      diffDiffRhoAtGaussPts(diff_diff_rho_at_gauss_pts),
      singularDispl(singular_displ), tetsInBlock(tets_in_block), rHo0(rho_0),
      boneN(bone_n) {
  if (forces_only_on_entities_row)
    forcesOnlyOnEntitiesRow = *forces_only_on_entities_row;
}
 
MoFEMErrorCode OpLhsBoneExplicitDerivariveWithHooke_dX::doWork(
    int row_side, int col_side, EntityType row_type, EntityType col_type,
    DataForcesAndSourcesCore::EntData &row_data,
    DataForcesAndSourcesCore::EntData &col_data) {
  MoFEMFunctionBegin;
  if (tetsInBlock.find(getFEEntityHandle()) == tetsInBlock.end()) {
    MoFEMFunctionReturnHot(0);
  }
 
  const int row_nb_dofs = row_data.getIndices().size();
  if (!row_nb_dofs)
    MoFEMFunctionReturnHot(0);
 
  const int col_nb_dofs = col_data.getIndices().size();
  if (!col_nb_dofs)
    MoFEMFunctionReturnHot(0);
 
  nA.resize(row_nb_dofs, col_nb_dofs, false);
  nA.clear();
 
  const int nb_gauss_pts = row_data.getN().size1();
  auto t_invH = getFTensor2FromMat<3, 3>(*commonData.invHMat);
  auto t_h = getFTensor2FromMat<3, 3>(*commonData.hMat);
  auto t_H = getFTensor2FromMat<3, 3>(*commonData.HMat);
  auto rho = getFTensor0FromVec(*rhoAtGaussPts);
  auto t_diff_rho = getFTensor1FromMat<3>(*diffRhoAtGaussPts);
  auto t_diff_diff_rho = getFTensor2FromMat<3, 3>(*diffDiffRhoAtGaussPts);
 
  auto t_cauchy_stress =
      getFTensor2SymmetricFromMat<3>(*commonData.cauchyStressMat);
 
  auto energy = getFTensor0FromVec(*commonData.energyVec);
  auto det_H = getFTensor0FromVec(*commonData.detHVec);
 
  auto t_initial_singular_displacement = getFTensor1FromMat<3>(singularDispl);
  auto t_row_diff_n = row_data.getFTensor1DiffN<3>();
 
  const int nb_row_base = row_data.getN().size2();
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
  FTensor::Index<'m', 3> m;
  FTensor::Index<'n', 3> n;
 
  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
    const double val = getVolume() * getGaussPts()(3, gg) * det_H;
 
    const double drho_coeff = (boneN / rho);
    const double denergy_drho = energy * drho_coeff;
    const double denergy_drho_drho =
        energy * (boneN - 1.) * boneN / (rho * rho);
    FTensor::Tensor4<double, 3, 3, 3, 3> t_F_dX;
    t_F_dX(i, j, k, l) = -(t_h(i, m) * t_invH(m, k)) * t_H(l, j);
 
    FTensor::Tensor2<double, 3, 3> t_energy_dX;
    t_energy_dX(k, l) = t_F_dX(i, j, k, l) * t_cauchy_stress(i, j);
 
    FTensor::Tensor1<double, 3> t_b;
    FTensor::Tensor1<double, 3> t_rho_dX_pulled;
    FTensor::Tensor1<double, 3> t_psi_dX_rho_dX;
    FTensor::Tensor1<double, 3> t_psi_dX_rho_dX_pulled;
    FTensor::Tensor2<double, 3, 3> t_rho_dX_dX_pulled;
    FTensor::Tensor2<double, 3, 3> t_energy_dX_pulled;
    t_rho_dX_pulled(i) = t_diff_rho(j) * t_invH(j, i);
    t_psi_dX_rho_dX(i) =
        t_energy_dX(i, j) * t_rho_dX_pulled(j) * val * drho_coeff;
    t_psi_dX_rho_dX_pulled(j) = t_psi_dX_rho_dX(i) * t_invH(i, j);
 
    auto t_row_base = row_data.getFTensor0N(gg, 0);
    int rr = 0;
    for (; rr != row_nb_dofs / 3; rr++) {
      FTensor::Tensor2<double *, 3, 3> t_mat(
          &nA(3 * rr + 0, 0), &nA(3 * rr + 0, 1), &nA(3 * rr + 0, 2),
          &nA(3 * rr + 1, 0), &nA(3 * rr + 1, 1), &nA(3 * rr + 1, 2),
          &nA(3 * rr + 2, 0), &nA(3 * rr + 2, 1), &nA(3 * rr + 2, 2), 3);
 
      FTensor::Tensor2<double, 3, 3> t_a;
      t_a(i, j) = val * denergy_drho_drho * t_diff_rho(i) * t_diff_rho(j);
      t_a(i, j) += val * denergy_drho * t_diff_diff_rho(i, j);
      t_b(i) = t_row_base * val * denergy_drho * t_diff_rho(i);
 
      auto t_col_base = col_data.getFTensor0N(gg, 0);
      auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
 
      for (int cc = 0; cc != col_nb_dofs / 3; cc++) {
 
        t_mat(i, j) += t_row_base * t_a(i, j) * t_col_base;
        t_mat(i, j) += t_row_base * t_a(i, j) * t_col_diff_base(k) *
                       t_initial_singular_displacement(k);
 
        t_mat(i, j) += t_b(i) * (t_invH(l, j) * t_col_diff_base(l));
        t_mat(i, j) +=
            t_row_base * t_psi_dX_rho_dX_pulled(i) * t_col_diff_base(j);
        ++t_mat;
        ++t_col_base;
        ++t_col_diff_base;
      }
      ++t_row_base;
      ++t_row_diff_n;
    }
    for (; rr != nb_row_base; ++rr) {
      ++t_row_base;
      ++t_row_diff_n;
    }
 
    ++t_invH;
    ++t_H;
    ++rho;
    ++energy;
    ++det_H;
    ++t_diff_diff_rho;
    ++t_diff_rho;
    ++t_h;
    //  ++t_D;
    //  ++t_strain;
    ++t_cauchy_stress;
    ++t_initial_singular_displacement;
  }
 
  // Assembly
  int *indices_ptr = &row_data.getIndices()[0];
  if (!forcesOnlyOnEntitiesRow.empty()) {
    iNdices.resize(row_nb_dofs);
    noalias(iNdices) = row_data.getIndices();
    indices_ptr = &iNdices[0];
    auto &dofs = row_data.getFieldDofs();
    auto dit = dofs.begin();
    for (int ii = 0; dit != dofs.end(); dit++, ii++) {
      if (auto dof = (*dit)) {
        if (forcesOnlyOnEntitiesRow.find(dof->getEnt()) ==
            forcesOnlyOnEntitiesRow.end()) {
          iNdices[ii] = -1;
        }
      }
    }
  }
 
  CHKERR MatSetValues(getFEMethod()->snes_B, row_nb_dofs, indices_ptr,
                      col_nb_dofs, &*col_data.getIndices().begin(),
                      &*nA.data().begin(), ADD_VALUES);
  MoFEMFunctionReturn(0);
}
 
OpLhsBoneExplicitDerivariveWithHooke_dx::
    OpLhsBoneExplicitDerivariveWithHooke_dx(
        HookeElement::DataAtIntegrationPts &common_data,
        boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
        boost::shared_ptr<MatrixDouble> diff_rho_at_gauss_pts,
        boost::shared_ptr<MatrixDouble> diff_diff_rho_at_gauss_pts,
        MatrixDouble &singular_displ, Range tets_in_block, const double rho_0,
        const double bone_n, Range *forces_only_on_entities_row)
    : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
          "MESH_NODE_POSITIONS", "SPATIAL_POSITION", UserDataOperator::OPROWCOL,
          false),
      commonData(common_data), rhoAtGaussPts(rho_at_gauss_pts),
      diffRhoAtGaussPts(diff_rho_at_gauss_pts),
      diffDiffRhoAtGaussPts(diff_diff_rho_at_gauss_pts),
      singularDispl(singular_displ), tetsInBlock(tets_in_block), rHo0(rho_0),
      boneN(bone_n) {
  if (forces_only_on_entities_row)
    forcesOnlyOnEntitiesRow = *forces_only_on_entities_row;
}
 
MoFEMErrorCode OpLhsBoneExplicitDerivariveWithHooke_dx::doWork(
    int row_side, int col_side, EntityType row_type, EntityType col_type,
    DataForcesAndSourcesCore::EntData &row_data,
    DataForcesAndSourcesCore::EntData &col_data) {
  MoFEMFunctionBegin;
  if (tetsInBlock.find(getFEEntityHandle()) == tetsInBlock.end()) {
    MoFEMFunctionReturnHot(0);
  }
  const int row_nb_dofs = row_data.getIndices().size();
  if (!row_nb_dofs) {
    MoFEMFunctionReturnHot(0);
  }
  const int col_nb_dofs = col_data.getIndices().size();
  if (!col_nb_dofs) {
    MoFEMFunctionReturnHot(0);
  }
  nA.resize(row_nb_dofs, col_nb_dofs, false);
  nA.clear();
 
  const int nb_row_base = row_data.getN().size2();
  const int nb_gauss_pts = row_data.getN().size1();
 
  auto t_initial_singular_displacement = getFTensor1FromMat<3>(singularDispl);
  auto t_row_diff_n = row_data.getFTensor1DiffN<3>();
 
  auto t_H = getFTensor2FromMat<3, 3>(*(commonData.HMat));
  auto t_h = getFTensor2FromMat<3, 3>(*(commonData.hMat));
  auto det_H = getFTensor0FromVec(*commonData.detHVec);
  auto t_invH = getFTensor2FromMat<3, 3>(*commonData.invHMat);
  auto energy = getFTensor0FromVec(*commonData.energyVec);
 
  auto t_cauchy_stress =
      getFTensor2SymmetricFromMat<3>(*(commonData.cauchyStressMat));
  auto rho = getFTensor0FromVec(*rhoAtGaussPts);
 
  auto t_diff_rho = getFTensor1FromMat<3>(*diffRhoAtGaussPts);
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
  FTensor::Index<'m', 3> m;
  FTensor::Tensor1<double, 3> t_a;
 
  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
    double val = getVolume() * getGaussPts()(3, gg) * det_H;
 
    auto t_row_base = row_data.getFTensor0N(gg, 0);
    const double drho_coef = (boneN / rho);
 
    int rr = 0;
    for (; rr != row_nb_dofs / 3; ++rr) {
 
      auto t_col_n = col_data.getFTensor0N(gg, 0);
      auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
 
      FTensor::Tensor2<double *, 3, 3> t_mat(
          &nA(3 * rr + 0, 0), &nA(3 * rr + 0, 1), &nA(3 * rr + 0, 2),
          &nA(3 * rr + 1, 0), &nA(3 * rr + 1, 1), &nA(3 * rr + 1, 2),
          &nA(3 * rr + 2, 0), &nA(3 * rr + 2, 1), &nA(3 * rr + 2, 2), 3);
      FTensor::Tensor2<double, 3, 3> t_b;
 
      t_a(i) =
          t_row_base * t_cauchy_stress(i, j) * t_diff_rho(j) * val * drho_coef;
 
      for (int cc = 0; cc != col_nb_dofs / 3; cc++) {
 
        t_mat(i, j) += t_a(i) * (t_invH(l, j) * t_col_diff_base(l));
 
        ++t_mat;
        ++t_col_diff_base;
        ++t_col_n;
      }
      ++t_row_base;
      ++t_row_diff_n;
    }
 
    for (; rr != nb_row_base; ++rr) {
      ++t_row_base;
      ++t_row_diff_n;
    }
 
    ++t_h;
    ++t_H;
    ++t_invH;
    ++rho;
    ++t_diff_rho;
    ++t_cauchy_stress;
    ++det_H;
    ++t_initial_singular_displacement;
  }
 
  int *indices_ptr = &row_data.getIndices()[0];
  if (!forcesOnlyOnEntitiesRow.empty()) {
    iNdices.resize(row_nb_dofs);
    noalias(iNdices) = row_data.getIndices();
    indices_ptr = &iNdices[0];
    auto &dofs = row_data.getFieldDofs();
    auto dit = dofs.begin();
    for (int ii = 0; dit != dofs.end(); dit++, ii++) {
      if(auto dof = (*dit)) {
        if (forcesOnlyOnEntitiesRow.find(dof->getEnt()) ==
            forcesOnlyOnEntitiesRow.end()) {
          iNdices[ii] = -1;
        }
      }
    }
  }
 
  CHKERR MatSetValues(getFEMethod()->snes_B, row_nb_dofs, indices_ptr,
                      col_nb_dofs, &*col_data.getIndices().begin(),
                      &*nA.data().begin(), ADD_VALUES);
 
  MoFEMFunctionReturn(0);
}
 
OpAnalyticalSpatialTraction::OpAnalyticalSpatialTraction(
    MoFEM::Interface &m_field, const bool &set_singular_coordinates,
    const Range &crack_front_nodes, const Range &crack_front_nodes_edges,
    const Range &crack_front_elements, PetscBool add_singularity,
    const std::string field_name, const Range tris,
    boost::ptr_vector<MethodForForceScaling> &methods_op,
    boost::ptr_vector<NeumannForcesSurface::MethodForAnalyticalForce>
        &analytical_force_op)
    : FaceElementForcesAndSourcesCore::UserDataOperator(
          field_name, UserDataOperator::OPROW),
      tRis(tris), methodsOp(methods_op), analyticalForceOp(analytical_force_op),
      volSideFe(m_field, set_singular_coordinates, crack_front_nodes,
                crack_front_nodes_edges, crack_front_elements, add_singularity),
      setSingularCoordinates(set_singular_coordinates) {}
 
MoFEMErrorCode
OpAnalyticalSpatialTraction::doWork(int side, EntityType type,
                                    DataForcesAndSourcesCore::EntData &data) {
  MoFEMFunctionBeginHot;
 
  if (data.getIndices().size() == 0)
    MoFEMFunctionReturnHot(0);
  EntityHandle ent = getNumeredEntFiniteElementPtr()->getEnt();
  if (tRis.find(ent) == tRis.end())
    MoFEMFunctionReturnHot(0);
 
  if (type == MBVERTEX) {
    CHKERR loopSideVolumes("ELASTIC", volSideFe);
  }
 
  int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
  int nb_row_dofs = data.getIndices().size() / rank;
 
  Nf.resize(data.getIndices().size(), false);
  Nf.clear();
 
  VectorDouble3 coords(3);
  VectorDouble3 normal(3);
  VectorDouble3 force(3);
  VectorDouble3 FTforce(3);
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
 
  //  If singular element apply transformation to h tensor
  if (type == MBVERTEX && volSideFe.singularElement) {
    double *inv_jac_ptr = &*volSideFe.invSJac.data().begin();
    FTensor::Tensor2<double *, 3, 3> t_inv_s_jac(
        inv_jac_ptr, &inv_jac_ptr[1], &inv_jac_ptr[2], &inv_jac_ptr[3],
        &inv_jac_ptr[4], &inv_jac_ptr[5], &inv_jac_ptr[6], &inv_jac_ptr[7],
        &inv_jac_ptr[8], 9);
    FTensor::Tensor1<double *, 3> t_normal = getFTensor1Normal();
    double nrm = 2 * getArea();
    FTensor::Tensor1<double, 3> t_s_normal;
    sNrm.resize(data.getN().size1(), false);
    for (unsigned int gg = 0; gg != data.getN().size1(); gg++) {
      // Push normal to tranformed singular element
      double det = volSideFe.detS[gg];
      t_s_normal(i) = det * t_inv_s_jac(j, i) * t_normal(j) / nrm;
      sNrm[gg] = sqrt(t_s_normal(i) * t_s_normal(i));
      ++t_inv_s_jac;
    }
  }
 
  FTensor::Tensor1<double *, 3> t_force(&force[0], &force[1], &force[2]);
 
  const auto &coords_at_gauss_pts = getCoordsAtGaussPts();
  const auto &normals_at_gauss_pts = getNormalsAtGaussPts();
 
  for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
 
    // get integration weight and Jacobian of integration point (area of face)
    double val = getGaussPts()(2, gg);
    val *= 0.5 * cblas_dnrm2(3, &getNormalsAtGaussPts()(gg, 0), 1);
    for (int dd = 0; dd != 3; dd++) {
      coords[dd] = coords_at_gauss_pts(gg, dd);
      normal[dd] = normals_at_gauss_pts(gg, dd);
    }
    if (volSideFe.singularElement) {
      for (int dd = 0; dd != 3; dd++) {
        coords[dd] += volSideFe.singularDisp(gg, dd);
      }
      val *= sNrm[gg];
    }
 
    for (boost::ptr_vector<
             NeumannForcesSurface::MethodForAnalyticalForce>::iterator vit =
             analyticalForceOp.begin();
         vit != analyticalForceOp.end(); vit++) {
 
      // Calculate force
      CHKERR vit->getForce(ent, coords, normal, force);
 
      for (int rr = 0; rr != 3; rr++) {
        cblas_daxpy(nb_row_dofs, val * force[rr], &data.getN()(gg, 0), 1,
                    &Nf[rr], rank);
      }
    }
  }
 
  // Scale force using user defined scaling operator
  CHKERR MethodForForceScaling::applyScale(getFEMethod(), methodsOp, Nf);
 
  // Assemble force into vector
  CHKERR VecSetValues(getFEMethod()->snes_f, data.getIndices().size(),
                      &data.getIndices()[0], &Nf[0], ADD_VALUES);
 
  MoFEMFunctionReturnHot(0);
}
 
OpAnalyticalMaterialTraction::OpAnalyticalMaterialTraction(
    MoFEM::Interface &m_field, const bool &set_singular_coordinates,
    const Range &crack_front_nodes, const Range &crack_front_nodes_edges,
    const Range &crack_front_elements, PetscBool add_singularity,
    const std::string field_name, const Range tris,
    boost::ptr_vector<MethodForForceScaling> &methods_op,
    boost::ptr_vector<NeumannForcesSurface::MethodForAnalyticalForce>
        &analytical_force_op,
    Range *forces_only_on_entities_row)
    : FaceElementForcesAndSourcesCore::UserDataOperator(
          field_name, UserDataOperator::OPROW),
      tRis(tris), methodsOp(methods_op), analyticalForceOp(analytical_force_op),
      volSideFe(m_field, set_singular_coordinates, crack_front_nodes,
                crack_front_nodes_edges, crack_front_elements, add_singularity),
      setSingularCoordinates(set_singular_coordinates) {
  hG = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
  HG = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
  volSideFe.getOpPtrVector().push_back(
      new OpCalculateVectorFieldGradient<3, 3>("SPATIAL_POSITION", hG));
  volSideFe.getOpPtrVector().push_back(
      new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS", HG));
  volSideFe.meshPositionsFieldName = "NONE";
  if (forces_only_on_entities_row) {
    forcesOnlyOnEntitiesRow = *forces_only_on_entities_row;
  }
}
 
MoFEMErrorCode
OpAnalyticalMaterialTraction::doWork(int side, EntityType type,
                                     DataForcesAndSourcesCore::EntData &data) {
  MoFEMFunctionBeginHot;
 
  if (data.getIndices().size() == 0)
    MoFEMFunctionReturnHot(0);
  EntityHandle ent = getNumeredEntFiniteElementPtr()->getEnt();
  if (tRis.find(ent) == tRis.end())
    MoFEMFunctionReturnHot(0);
 
  if (type == MBVERTEX) {
    CHKERR loopSideVolumes("MATERIAL", volSideFe);
  }
  if (type != MBVERTEX) {
    if (volSideFe.singularElement) {
      const EntityHandle ent = data.getFieldDofs()[0]->getEnt();
      int side_number, sense, offset;
      CHKERR volSideFe.mField.get_moab().side_number(
          volSideFe.numeredEntFiniteElementPtr->getEnt(), ent, side_number,
          sense, offset);
      data.getN() =
          volSideFe.getEntData(H1, type, side_number).getN(data.getBase());
    }
  } else if (volSideFe.singularElement) {
    // Set value of base functions at singular integration points
    // FIXME: That need to be moved to separate operator such that values
    // of base function are appropriately evaluated
    for (unsigned int gg = 0; gg != volSideFe.gaussPts.size2(); gg++) {
      for (int dd = 0; dd != 3; dd++) {
        volSideFe.gaussPts(dd, gg) += volSideFe.singularRefCoords(gg, dd);
      }
    }
    CHKERR ShapeMBTET(
        &*volSideFe.getEntData(H1, MBVERTEX, 0).getN(NOBASE).data().begin(),
        &volSideFe.gaussPts(0, 0), &volSideFe.gaussPts(1, 0),
        &volSideFe.gaussPts(2, 0), volSideFe.gaussPts.size2());
    CHKERR volSideFe.calHierarchicalBaseFunctionsOnElement(data.getBase());
    map<int, int> map_face_nodes;
    for (unsigned int gg = 0; gg != data.getN().size1(); gg++) {
      for (int nn = 0; nn != 3; nn++) {
        data.getN()(gg, nn) =
            volSideFe.getEntData(H1, MBVERTEX, 0)
                .getN(data.getBase())(gg, volSideFe.getFaceConnMap()[nn]);
      }
    }
  }
 
  int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
  int nb_row_base = data.getIndices().size() / rank;
 
  Nf.resize(data.getIndices().size(), false);
  Nf.clear();
 
  VectorDouble3 coords(3);
  VectorDouble3 normal(3);
  VectorDouble3 force(3);
  VectorDouble3 FTforce(3);
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
 
  //  If singular element apply transformation to h tensor
  if (type == MBVERTEX && volSideFe.singularElement) {
    //  cerr << volSideFe.sJac << endl;
    //  cerr << volSideFe.invSJac << endl;
    double *inv_jac_ptr = &*volSideFe.invSJac.data().begin();
    FTensor::Tensor2<double *, 3, 3> t_inv_s_jac(
        inv_jac_ptr, &inv_jac_ptr[1], &inv_jac_ptr[2], &inv_jac_ptr[3],
        &inv_jac_ptr[4], &inv_jac_ptr[5], &inv_jac_ptr[6], &inv_jac_ptr[7],
        &inv_jac_ptr[8], 9);
    auto t_hg = getFTensor2FromMat<3, 3>(*hG);
    auto t_normal = getFTensor1Normal();
    double nrm = 2 * getArea();
    FTensor::Tensor1<double, 3> t_s_normal;
    sNrm.resize(data.getN().size1(), false);
    for (unsigned int gg = 0; gg != data.getN().size1(); gg++) {
      // Push normal to transformed singular element
      double det = volSideFe.detS[gg];
      t_s_normal(i) = det * t_inv_s_jac(j, i) * t_normal(j) / nrm;
      sNrm[gg] = sqrt(t_s_normal(i) * t_s_normal(i));
      FTensor::Tensor2<double, 3, 3> t_tmp;
      t_tmp(i, j) = t_hg(i, k) * t_inv_s_jac(k, j);
      t_hg(i, j) = t_tmp(i, j);
      ++t_inv_s_jac;
      ++t_hg;
    }
  }
 
  FTensor::Tensor1<double *, 3> t_force(&force[0], &force[1], &force[2]);
  FTensor::Tensor1<double *, 3> t_ft_force(&FTforce[0], &FTforce[1],
                                           &FTforce[2]);
  auto t_hg = getFTensor2FromMat<3, 3>(*hG);
  auto t_Hg = getFTensor2FromMat<3, 3>(*HG);
  FTensor::Tensor2<double, 3, 3> t_inv_Hg;
  FTensor::Tensor2<double, 3, 3> t_F;
 
  //  cerr << "\n\n";
  for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
 
    // get integration weight and Jacobian of integration point (area of face)
    double val = getGaussPts()(2, gg);
    val *= 0.5 * cblas_dnrm2(3, &getNormalsAtGaussPts()(gg, 0), 1);
    for (int dd = 0; dd != 3; dd++) {
      coords[dd] = getCoordsAtGaussPts()(gg, dd);
      normal[dd] = getNormalsAtGaussPts()(gg, dd);
    }
    if (volSideFe.singularElement) {
      for (int dd = 0; dd != 3; dd++) {
        coords[dd] += volSideFe.singularDisp(gg, dd);
      }
      val *= sNrm[gg];
    }
 
    for (boost::ptr_vector<
             NeumannForcesSurface::MethodForAnalyticalForce>::iterator vit =
             analyticalForceOp.begin();
         vit != analyticalForceOp.end(); vit++) {
 
      // Calculate force
      CHKERR vit->getForce(ent, coords, normal, force);
      // Invert matrix H
      double det;
      CHKERR determinantTensor3by3(t_Hg, det);
      CHKERR invertTensor3by3(t_Hg, det, t_inv_Hg);
 
      // Calculate gradient of deformation
      t_F(i, k) = t_hg(i, j) * t_inv_Hg(j, k);
      t_ft_force(i) = -t_F(j, i) * t_force(j);
 
      for (int rr = 0; rr != 3; rr++) {
        cblas_daxpy(nb_row_base, val * FTforce[rr], &data.getN()(gg, 0), 1,
                    &Nf[rr], rank);
      }
    }
 
    ++t_hg;
    ++t_Hg;
  }
 
  // Scale force using user defined scaling operator
  CHKERR MethodForForceScaling::applyScale(getFEMethod(), methodsOp, Nf);
 
  // Assemble force into vector
  int nb_dofs = data.getIndices().size();
  int *indices_ptr = &data.getIndices()[0];
  if (!forcesOnlyOnEntitiesRow.empty()) {
    iNdices.resize(nb_dofs, false);
    noalias(iNdices) = data.getIndices();
    indices_ptr = &iNdices[0];
    auto &dofs = data.getFieldDofs();
    auto dit = dofs.begin();
    for (int ii = 0; dit != dofs.end(); dit++, ii++) {
      if (auto dof = (*dit)) {
        if (forcesOnlyOnEntitiesRow.find(dof->getEnt()) ==
            forcesOnlyOnEntitiesRow.end()) {
          iNdices[ii] = -1;
        }
      }
    }
  }
  CHKERR VecSetValues(getFEMethod()->snes_f, nb_dofs, indices_ptr, &Nf[0],
                      ADD_VALUES);
  MoFEMFunctionReturnHot(0);
}
 
OpAleLhsWithDensitySingularElement_dX_dX::
    OpAleLhsWithDensitySingularElement_dX_dX(
        const std::string row_field, const std::string col_field,
        boost::shared_ptr<HookeElement::DataAtIntegrationPts> &data_at_pts,
        boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
        boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts,
        const double rho_n, const double rho_0,
        boost::shared_ptr<MatrixDouble> singular_displacement)
    : UserDataOperator(row_field, col_field, OPROWCOL, true),
      dataAtPts(data_at_pts), rhoAtGaussPtsPtr(rho_at_gauss_pts),
      rhoGradAtGaussPtsPtr(rho_grad_at_gauss_pts),
      singularDisplacement(singular_displacement), rhoN(rho_n), rHo0(rho_0) {}
 
MoFEMErrorCode OpAleLhsWithDensitySingularElement_dX_dX::doWork(
    int row_side, int col_side, EntityType row_type, EntityType col_type,
    DataForcesAndSourcesCore::EntData &row_data,
    DataForcesAndSourcesCore::EntData &col_data) {
 
  MoFEMFunctionBegin;
  nbRows = row_data.getIndices().size();
  if (!nbRows)
    MoFEMFunctionReturnHot(0);
 
  nbCols = col_data.getIndices().size();
  if (!nbCols)
    MoFEMFunctionReturnHot(0);
 
  K.resize(nbRows, nbCols, false);
  K.clear();
 
  nbIntegrationPts = getGaussPts().size2();
  if (row_side == col_side && row_type == col_type) {
    isDiag = true;
  } else {
    isDiag = false;
  }
 
  CHKERR iNtegrate(row_data, col_data);
  CHKERR aSsemble(row_data, col_data);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpAleLhsWithDensitySingularElement_dX_dX::iNtegrate(
    DataForcesAndSourcesCore::EntData &row_data,
    DataForcesAndSourcesCore::EntData &col_data) {
  MoFEMFunctionBegin;
  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
    return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(
        &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
        &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
        &m(r + 2, c + 2));
  };
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
 
  double vol = getVolume();
 
  auto t_w = getFTensor0IntegrationWeight();
 
  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
  const int row_nb_base_fun = row_data.getN().size2();
 
  auto rho = getFTensor0FromVec(*rhoAtGaussPtsPtr);
  auto t_grad_rho = getFTensor1FromMat<3>(*rhoGradAtGaussPtsPtr);
 
  auto t_eshelby_stress =
      getFTensor2FromMat<3, 3>(*dataAtPts->eshelbyStressMat);
  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
  auto &det_H = *dataAtPts->detHVec;
 
  auto t_initial_singular_displacement =
      getFTensor1FromMat<3>(*singularDisplacement);
 
  for (int gg = 0; gg != nbIntegrationPts; ++gg) {
 
    double a = t_w * vol * det_H[gg];
    const double stress_dho_coef = (rhoN / rho);
 
    int rr = 0;
    for (; rr != nbRows / 3; ++rr) {
 
      auto t_m = get_tensor2(K, 3 * rr, 0);
 
      FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
      t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
 
      FTensor::Tensor1<double, 3> t_row_stress;
      t_row_stress(i) = a * t_row_diff_base_pulled(j) * t_eshelby_stress(i, j);
      FTensor::Tensor2<double, 3, 3> t_a;
      t_a(i, k) = t_row_stress(i) * stress_dho_coef * t_grad_rho(k);
 
      auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
      auto t_col_base = col_data.getFTensor0N(gg, 0);
      for (int cc = 0; cc != nbCols / 3; ++cc) {
 
        t_m(i, k) +=
            t_a(i, k) * t_col_diff_base(l) * t_initial_singular_displacement(l);
 
        ++t_col_base;
        ++t_col_diff_base;
        ++t_m;
      }
      ++t_row_diff_base;
    }
 
    for (; rr != row_nb_base_fun; ++rr)
      ++t_row_diff_base;
 
    ++t_initial_singular_displacement;
    ++t_w;
    ++t_eshelby_stress;
    ++t_invH;
    ++rho;
    ++t_grad_rho;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpAleLhsWithDensitySingularElement_dX_dX::aSsemble(
    DataForcesAndSourcesCore::EntData &row_data,
    DataForcesAndSourcesCore::EntData &col_data) {
  MoFEMFunctionBegin;
 
  const int *row_indices = &*row_data.getIndices().data().begin();
  const int *col_indices = &*col_data.getIndices().data().begin();
 
  auto &data = *dataAtPts;
  Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
                                             : getFEMethod()->snes_B;
  // assemble local matrix
  CHKERR MatSetValues(B, nbRows, row_indices, nbCols, col_indices,
                      &*K.data().begin(), ADD_VALUES);
 
  if (!isDiag && sYmm) {
    // if not diagonal term and since global matrix is symmetric assemble
    // transpose term.
    transK.resize(K.size2(), K.size1(), false);
    noalias(transK) = trans(K);
    CHKERR MatSetValues(B, nbCols, col_indices, nbRows, row_indices,
                        &*transK.data().begin(), ADD_VALUES);
  }
  MoFEMFunctionReturn(0);
}
 
OpAleLhsWithDensitySingularElement_dx_dX::
    OpAleLhsWithDensitySingularElement_dx_dX(
        const std::string row_field, const std::string col_field,
        boost::shared_ptr<HookeElement::DataAtIntegrationPts> &data_at_pts,
        boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
        boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts,
        const double rho_n, const double rho_0,
        boost::shared_ptr<MatrixDouble> singular_displacement)
    : UserDataOperator(row_field, col_field, OPROWCOL, false),
      dataAtPts(data_at_pts), rhoAtGaussPtsPtr(rho_at_gauss_pts),
      rhoGradAtGaussPtsPtr(rho_grad_at_gauss_pts),
      singularDisplacement(singular_displacement), rhoN(rho_n), rHo0(rho_0) {}
 
MoFEMErrorCode OpAleLhsWithDensitySingularElement_dx_dX::doWork(
    int row_side, int col_side, EntityType row_type, EntityType col_type,
    DataForcesAndSourcesCore::EntData &row_data,
    DataForcesAndSourcesCore::EntData &col_data) {
 
  MoFEMFunctionBegin;
  nbRows = row_data.getIndices().size();
  if (!nbRows)
    MoFEMFunctionReturnHot(0);
 
  nbCols = col_data.getIndices().size();
  if (!nbCols)
    MoFEMFunctionReturnHot(0);
 
  K.resize(nbRows, nbCols, false);
  K.clear();
 
  nbIntegrationPts = getGaussPts().size2();
  if (row_side == col_side && row_type == col_type) {
    isDiag = true;
  } else {
    isDiag = false;
  }
 
  CHKERR iNtegrate(row_data, col_data);
  CHKERR aSsemble(row_data, col_data);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpAleLhsWithDensitySingularElement_dx_dX::iNtegrate(
    DataForcesAndSourcesCore::EntData &row_data,
    DataForcesAndSourcesCore::EntData &col_data) {
  MoFEMFunctionBegin;
  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
    return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(
        &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
        &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
        &m(r + 2, c + 2));
  };
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
 
  double vol = getVolume();
 
  auto t_w = getFTensor0IntegrationWeight();
 
  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
  const int row_nb_base_fun = row_data.getN().size2();
 
  auto rho = getFTensor0FromVec(*rhoAtGaussPtsPtr);
  auto t_grad_rho = getFTensor1FromMat<3>(*rhoGradAtGaussPtsPtr);
 
  auto t_cauchy_stress =
      getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
  auto &det_H = *dataAtPts->detHVec;
 
  auto t_initial_singular_displacement =
      getFTensor1FromMat<3>(*singularDisplacement);
 
  for (int gg = 0; gg != nbIntegrationPts; ++gg) {
 
    double a = t_w * vol * det_H[gg];
    const double stress_dho_coef = (rhoN / rho);
 
    int rr = 0;
    for (; rr != nbRows / 3; ++rr) {
 
      auto t_m = get_tensor2(K, 3 * rr, 0);
 
      FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
      t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
 
      FTensor::Tensor1<double, 3> t_row_stress;
      t_row_stress(i) = a * t_row_diff_base_pulled(j) * t_cauchy_stress(i, j);
      FTensor::Tensor2<double, 3, 3> t_a;
      t_a(i, k) = t_row_stress(i) * stress_dho_coef * t_grad_rho(k);
 
      auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
      auto t_col_base = col_data.getFTensor0N(gg, 0);
      for (int cc = 0; cc != nbCols / 3; ++cc) {
 
        t_m(i, k) +=
            t_a(i, k) * t_col_diff_base(l) * t_initial_singular_displacement(l);
 
        ++t_col_base;
        ++t_col_diff_base;
        ++t_m;
      }
      ++t_row_diff_base;
    }
 
    for (; rr != row_nb_base_fun; ++rr)
      ++t_row_diff_base;
 
    ++t_initial_singular_displacement;
    ++t_w;
    ++t_cauchy_stress;
    ++t_invH;
    ++rho;
    ++t_grad_rho;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpAleLhsWithDensitySingularElement_dx_dX::aSsemble(
    DataForcesAndSourcesCore::EntData &row_data,
    DataForcesAndSourcesCore::EntData &col_data) {
  MoFEMFunctionBegin;
 
  const int *row_indices = &*row_data.getIndices().data().begin();
  const int *col_indices = &*col_data.getIndices().data().begin();
 
  auto &data = *dataAtPts;
  Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
                                             : getFEMethod()->snes_B;
  // assemble local matrix
  CHKERR MatSetValues(B, nbRows, row_indices, nbCols, col_indices,
                      &*K.data().begin(), ADD_VALUES);
 
  if (!isDiag && sYmm) {
    // if not diagonal term and since global matrix is symmetric assemble
    // transpose term.
    transK.resize(K.size2(), K.size1(), false);
    noalias(transK) = trans(K);
    CHKERR MatSetValues(B, nbCols, col_indices, nbRows, row_indices,
                        &*transK.data().begin(), ADD_VALUES);
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
OpGetDensityFieldForTesting::doWork(int row_side, EntityType row_type,
                                    HookeElement::EntData &row_data) {
  MoFEMFunctionBegin;
  if (row_type != MBVERTEX)
    MoFEMFunctionReturnHot(0);
  // get number of integration points
  const int nb_integration_pts = getGaussPts().size2();
  rhoAtGaussPtsPtr->resize(nb_integration_pts, false);
  rhoAtGaussPtsPtr->clear();
  rhoGradAtGaussPtsPtr->resize(3, nb_integration_pts, false);
  rhoGradAtGaussPtsPtr->clear();
  rhoGradGradAtGaussPtsPtr->resize(9, nb_integration_pts, false);
  rhoGradGradAtGaussPtsPtr->clear();
 
  auto t_rho = getFTensor0FromVec(*rhoAtGaussPtsPtr);
  auto t_grad_rho = getFTensor1FromMat<3>(*rhoGradAtGaussPtsPtr);
  auto t_grad_grad_rho = getFTensor2FromMat<3, 3>(*rhoGradGradAtGaussPtsPtr);
 
  MatrixDouble &coords = *matCoordsPtr;
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
 
  singularDisp.resize(3, nb_integration_pts, false);
  singularDisp.clear();
 
  auto t_material_positions = getFTensor1FromMat<3>(coords);
  auto t_mwls_material_positions = getFTensor1FromMat<3>(coords);
  if (feSingularPtr->singularElement) {
    MatrixDouble singular_current_displacement;
 
    singular_current_displacement.resize(3, nb_integration_pts, false);
    singular_current_displacement.clear();
 
    FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_singular_displacement(
        &feSingularPtr->singularDisp(0, 0), &feSingularPtr->singularDisp(0, 1),
        &feSingularPtr->singularDisp(0, 2));
    auto t_inital_singular_displacement = getFTensor1FromMat<3>(singularDisp);
    auto t_current_singular_displacement =
        getFTensor1FromMat<3>(singular_current_displacement);
    if (!matGradPosAtPtsPtr)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Matrix for gradient of positions not allocated");
    auto t_H = getFTensor2FromMat<3, 3>(*matGradPosAtPtsPtr);
    for (int gg = 0; gg != nb_integration_pts; gg++) {
 
      t_inital_singular_displacement(i) = t_singular_displacement(i);
      t_current_singular_displacement(i) =
          t_H(i, j) * t_singular_displacement(j);
 
      t_material_positions(i) += t_current_singular_displacement(i);
 
      ++t_material_positions;
      ++t_H;
      ++t_singular_displacement;
      ++t_inital_singular_displacement;
      ++t_current_singular_displacement;
    }
  }
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    t_rho =
        1 + t_mwls_material_positions(i) *
                t_mwls_material_positions(i); // density is equal to
                                              // 1+x^2+y^2+z^2 (x,y,z coordines)
    t_grad_rho(i) = 2 * t_mwls_material_positions(i);
    t_grad_grad_rho(0, 0) = 2.;
    t_grad_grad_rho(1, 1) = 2.;
    t_grad_grad_rho(2, 2) = 2.;
 
    ++t_rho;
    ++t_mwls_material_positions;
    ++t_grad_rho;
    ++t_grad_grad_rho;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
OpSetTagRangeOnSkin::doWork(int side, EntityType type,
                               DataForcesAndSourcesCore::EntData &data) {
  MoFEMFunctionBeginHot;
 
  if (data.getFieldData().size() == 0)
    MoFEMFunctionReturnHot(0);
 
  int tag_val = tagValue;
  EntityHandle tri_ent = getNumeredEntFiniteElementPtr()->getEnt();
  if (rangeToTag.find(tri_ent) == rangeToTag.end())
    tag_val = 0;
 
  int def_VAL = 0;
  Tag th;
  CHKERR postProcMesh.tag_get_handle(tagName.c_str(), 1, MB_TYPE_INTEGER, th,
                                     MB_TAG_CREAT | MB_TAG_SPARSE, &def_VAL);
 
  const int nb_gauss_pts = mapGaussPts.size();
  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
 
    // check for existing mark
    int check_data;
    CHKERR postProcMesh.tag_get_data(th, &mapGaussPts[gg], 1, &check_data);
 
    if (check_data == 0)
      CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, &tag_val);
  }
  MoFEMFunctionReturnHot(0);
}
 
} // namespace FractureMechanics