CrackFrontElement.hpp

Go to the documentation of this file.

/** \file CrackFrontElement.hpp
 */
 
/* 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/>. */
 
#ifndef __CRACK_FRONT_ELEMENT_HPP__
#define __CRACK_FRONT_ELEMENT_HPP__
 
namespace FractureMechanics {
 
enum Materials { HOOKE, KIRCHHOFF, NEOHOOKEAN, BONEHOOKE, LASTOP };
 
template <typename E, typename B> class TwoType {};
 
/**
 *
 * Use idea from \cite barsoum1976use to approximate singularity at crack tip.
 * The idea is to shift mid noted to quarter edge length to crack tip on edges
 * adjacent to crack front.
 *
 */
template <typename ELEMENT, typename BASE_ELEMENT>
struct CrackFrontSingularBase : public ELEMENT {
 
  Range crackFrontNodes;
  Range crackFrontNodesEdges;
  Range crackFrontElements;
 
  MatrixDouble sJac;
  MatrixDouble invSJac;
  MatrixDouble singularDisp;
  MatrixDouble singularRefCoords;
  VectorDouble detS;
 
  bool singularElement;
  bool singularEdges[6];
  bool singularNodes[4];
 
  MoFEMErrorCode getElementOptions();
 
  CrackFrontSingularBase(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)
      : ELEMENT(m_field), crackFrontNodes(crack_front_nodes),
        crackFrontNodesEdges(crack_front_nodes_edges),
        crackFrontElements(crack_front_elements),
        addSingularity(add_singularity),
        setSingularCoordinatesPriv(set_singular_coordinates) {
    ierr = getElementOptions();
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
  }
 
  bool setSingularCoordinates;
  CrackFrontSingularBase(MoFEM::Interface &m_field)
      : ELEMENT(m_field), setSingularCoordinatesPriv(setSingularCoordinates) {}
 
  virtual ~CrackFrontSingularBase() = default;
 
  MoFEMErrorCode operator()() {
    MoFEMFunctionBegin;
 
    if (this->numeredEntFiniteElementPtr->getEntType() != MBTET)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Element type not implemented");
 
    this->getElementPolynomialBase() =
        boost::shared_ptr<BaseFunction>(new TetPolynomialBase());
 
    CHKERR this->createDataOnElement(MBTET);
    CHKERR this->calculateVolumeAndJacobian();
    CHKERR this->getSpaceBaseAndOrderOnElement();
    CHKERR this->setIntegrationPts();
    CHKERR this->calculateCoordinatesAtGaussPts();
    CHKERR this->calHierarchicalBaseFunctionsOnElement();
    CHKERR this->calculateHOJacobian();
    CHKERR this->transformBaseFunctions();
    
    // Iterate over operators
    CHKERR this->loopOverOperators();
 
    MoFEMFunctionReturn(0);
  }
 
  PetscBool addSingularity;
  PetscBool refineIntegration;
  int refinementLevels;
 
  MoFEMErrorCode getSpaceBaseAndOrderOnElement() {
    return getSpaceBaseAndOrderOnElementImpl(TwoType<ELEMENT, BASE_ELEMENT>());
  }
 
  template <typename E, typename B>
  MoFEMErrorCode getSpaceBaseAndOrderOnElementImpl(TwoType<E, B>) {
    return E::getSpaceBaseAndOrderOnElement();
  }
 
  /// Partial specialization for volume elements
  template <typename E>
  MoFEMErrorCode getSpaceBaseAndOrderOnElementImpl(
      TwoType<E, VolumeElementForcesAndSourcesCore>) {
    MoFEMFunctionBeginHot;
    CHKERR E::getSpaceBaseAndOrderOnElement();
    // singularElement = false;
    // MoFEMFunctionReturnHot(0);
    // Determine if this element is singular
    EntityHandle ent = this->numeredEntFiniteElementPtr->getEnt();
    if (crackFrontElements.find(ent) != crackFrontElements.end()) {
      singularElement = true;
      for (int nn = 0; nn != 4; nn++) {
        if (crackFrontNodes.find(this->conn[nn]) != crackFrontNodes.end()) {
          singularNodes[nn] = true;
        } else {
          singularNodes[nn] = false;
        }
      }
      for (int ee = 0; ee != 6; ee++) {
        EntityHandle edge;
        CHKERR this->mField.get_moab().side_element(ent, 1, ee, edge);
        if (crackFrontNodesEdges.find(edge) != crackFrontNodesEdges.end()) {
          singularEdges[ee] = true;
        } else {
          singularEdges[ee] = false;
        }
      }
    } else {
      singularElement = false;
    }
    MoFEMFunctionReturnHot(0);
  }
 
  MoFEMErrorCode calculateHOJacobian() {
    return calculateHOJacobianImpl(
        TwoType<ELEMENT, VolumeElementForcesAndSourcesCore>());
  }
 
  /// Partial specialization for volume element
  template <typename E>
  MoFEMErrorCode
  calculateHOJacobianImpl(TwoType<E, VolumeElementForcesAndSourcesCore>) {
    MoFEMFunctionBegin;
 
    if (setSingularCoordinatesPriv) {
      singularElement = false;
      MoFEMFunctionReturnHot(0);
    }
    if (addSingularity == PETSC_FALSE) {
      singularElement = false;
      MoFEMFunctionReturnHot(0);
    }
 
    if (singularElement) {
 
      const int edge_nodes[6][2] = {{0, 1}, {1, 2}, {2, 0},
                                    {0, 3}, {1, 3}, {2, 3}};
      FTensor::Index<'i', 3> i;
      FTensor::Index<'j', 3> j;
      FTensor::Index<'k', 3> k;
 
      singularDisp.resize(this->coordsAtGaussPts.size1(),
                          this->coordsAtGaussPts.size2());
      singularDisp.clear();
      singularRefCoords.resize(this->coordsAtGaussPts.size1(),
                               this->coordsAtGaussPts.size2());
      singularRefCoords.clear();
 
      const size_t nb_gauss_pts = this->gaussPts.size2();
      sJac.resize(nb_gauss_pts, 9, false);
      double *s_jac_ptr = &sJac(0, 0);
      // Set jacobian to 1 on diagonal
      {
        FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_s_jac(
            s_jac_ptr, &s_jac_ptr[1], &s_jac_ptr[2], &s_jac_ptr[3],
            &s_jac_ptr[4], &s_jac_ptr[5], &s_jac_ptr[6], &s_jac_ptr[7],
            &s_jac_ptr[8]);
        for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
          t_s_jac(i, j) = 0;
          t_s_jac(0, 0) = 1;
          t_s_jac(1, 1) = 1;
          t_s_jac(2, 2) = 1;
          ++t_s_jac;
        }
      }
 
      // get base fuctions direvatives
      double diff_base_n[12];
      CHKERR ShapeDiffMBTET(diff_base_n);
      FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_diff_base_n(
          diff_base_n, &diff_base_n[1], &diff_base_n[2]);
      double diff_n[12];
      FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_diff_n(
          diff_n, &diff_n[1], &diff_n[2]);
      for (int nn = 0; nn != 4; nn++) {
        t_diff_n(i) = this->tInvJac(j, i) * t_diff_base_n(j);
        ++t_diff_n;
        ++t_diff_base_n;
      }
      MatrixDouble &shape_n =
          this->dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
 
      const double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
 
      // Calculate singular jacobian
      for (int ee = 0; ee != 6; ee++) {
        if (!singularEdges[ee])
          continue;
        // i0 is node at crack front, i1 is node inside body
        const int i0 = edge_nodes[ee][0];
        const int i1 = edge_nodes[ee][1];
        // Only one of the nodes at edge have to be at crack front
        if ((singularNodes[i0] && singularNodes[i1]) ||
            !(singularNodes[i0] || singularNodes[i1])) {
          PetscPrintf(PETSC_COMM_SELF, "Warning singular edge on both ends\n");
          // EntityHandle out_meshset;
          // CHKERR this->mField.get_moab().create_meshset(
          //     MESHSET_SET, out_meshset);
          // CHKERR this->mField.get_moab().add_entities(out_meshset,
          //                                             crackFrontNodesEdges);
          // CHKERR this->mField.get_moab().write_file(
          //     "debug_error_crack_front_nodes_edges.vtk", "VTK", "",
          //     &out_meshset, 1);
          // CHKERR this->mField.get_moab().delete_entities(&out_meshset, 1);
          // CHKERR this->mField.get_moab().create_meshset(
          //     MESHSET_SET, out_meshset);
          // CHKERR this->mField.get_moab().add_entities(out_meshset,
          //                                             crackFrontElements);
          // CHKERR this->mField.get_moab().write_file(
          //     "debug_error_crack_elements.vtk", "VTK", "", &out_meshset, 1);
          // CHKERR this->mField.get_moab().delete_entities(&out_meshset, 1);
          continue;
          // SETERRQ(
          //     PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
          //     "Huston we have problem, only one on the edge can be
          //     singular");
        }
        // edge directions oriented from crack front node towards inside
        FTensor::Tensor1<double, 3> t_base_dir(
            base_coords[3 * i1 + 0] - base_coords[3 * i0 + 0],
            base_coords[3 * i1 + 1] - base_coords[3 * i0 + 1],
            base_coords[3 * i1 + 2] - base_coords[3 * i0 + 2]);
        // set orientation
        if (singularNodes[edge_nodes[ee][0]]) {
          t_base_dir(i) *= -1;
        }
        // push direction to current configuration
        FTensor::Tensor1<double, 3> t_dir;
        t_dir(i) = this->tJac(i, j) * t_base_dir(j);
        FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_s_jac(
            s_jac_ptr, &s_jac_ptr[1], &s_jac_ptr[2], &s_jac_ptr[3],
            &s_jac_ptr[4], &s_jac_ptr[5], &s_jac_ptr[6], &s_jac_ptr[7],
            &s_jac_ptr[8]);
        FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>
            t_singular_displacement(&singularDisp(0, 0), &singularDisp(0, 1),
                                    &singularDisp(0, 2));
        FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>
            t_singular_ref_displacement(&singularRefCoords(0, 0),
                                        &singularRefCoords(0, 1),
                                        &singularRefCoords(0, 2));
        // Calculate jacobian
        for (size_t gg = 0; gg != nb_gauss_pts; gg++) {
          double t_n = shape_n(gg, i0) * shape_n(gg, i1);
          FTensor::Tensor1<double, 3> t_diff_n(
              shape_n(gg, i0) * diff_n[3 * i1 + 0] +
                  shape_n(gg, i1) * diff_n[3 * i0 + 0],
              shape_n(gg, i0) * diff_n[3 * i1 + 1] +
                  shape_n(gg, i1) * diff_n[3 * i0 + 1],
              shape_n(gg, i0) * diff_n[3 * i1 + 2] +
                  shape_n(gg, i1) * diff_n[3 * i0 + 2]);
          t_s_jac(i, j) += t_dir(i) * t_diff_n(j);
          t_singular_displacement(i) += t_dir(i) * t_n;
          t_singular_ref_displacement(i) += t_base_dir(i) * t_n;
          ++t_singular_displacement;
          ++t_singular_ref_displacement;
          ++t_s_jac;
        }
      }
 
      // invert singular jacobian and set integration points
      {
        FTensor::Tensor2<double *, 3, 3> t_s_jac(
            s_jac_ptr, &s_jac_ptr[1], &s_jac_ptr[2], &s_jac_ptr[3],
            &s_jac_ptr[4], &s_jac_ptr[5], &s_jac_ptr[6], &s_jac_ptr[7],
            &s_jac_ptr[8], 9);
        // Allocate memory for singular inverse jacobian and setup tensor.
        invSJac.resize(nb_gauss_pts, 9, false);
        double *inv_s_jac_ptr = &invSJac(0, 0);
        FTensor::Tensor2<double *, 3, 3> t_inv_s_jac(
            inv_s_jac_ptr, &inv_s_jac_ptr[1], &inv_s_jac_ptr[2],
            &inv_s_jac_ptr[3], &inv_s_jac_ptr[4], &inv_s_jac_ptr[5],
            &inv_s_jac_ptr[6], &inv_s_jac_ptr[7], &inv_s_jac_ptr[8], 9);
        // Sum of all determinist has to be 1
        // Calculate integrations weights
        detS.resize(nb_gauss_pts, false);
        for (size_t gg = 0; gg != nb_gauss_pts; gg++) {
          double det;
          CHKERR determinantTensor3by3(t_s_jac, det);
          CHKERR invertTensor3by3(t_s_jac, det, t_inv_s_jac);
          detS[gg] = det;
          ++t_inv_s_jac;
          ++t_s_jac;
        }
      }
    }
 
    MoFEMFunctionReturn(0);
  }
 
  /**
   * @return returning negative number -1, i.e. special user quadrature is
   * generated
   */
  int getRule(int order) {
    return getRuleImpl(TwoType<ELEMENT, BASE_ELEMENT>(), order);
  }
 
  // Generic specialization
  template <typename E, typename B> int getRuleImpl(TwoType<E, B>, int order) {
    return E::getRule(order);
  }
 
  /**
   * \brief Generate specific user integration quadrature
   *
   * Element is refined at the crack top and quadrature is set to each of
   * sub-tets.
   *
   */
  MoFEMErrorCode setGaussPts(int order) {
    return setGaussPtsImpl(TwoType<ELEMENT, BASE_ELEMENT>(), order);
  }
 
  // Generic specialization
  template <typename E, typename B>
  int setGaussPtsImpl(TwoType<E, B>, int order) {
    return E::setGaussPts(order);
  }
 
private:
  const bool &setSingularCoordinatesPriv;
 
  MatrixDouble refCoords;
  MatrixDouble refGaussCoords;
};
 
template <typename ELEMENT, typename BASE_ELEMENT>
MoFEMErrorCode
CrackFrontSingularBase<ELEMENT, BASE_ELEMENT>::getElementOptions() {
  MoFEMFunctionBeginHot;
  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Get singular element options",
                           "none");
  CHKERRQ(ierr);
  refinementLevels = 3;
  ierr = PetscOptionsInt("-se_number_of_refinement_levels",
                         "approximation geometry order", "", refinementLevels,
                         &refinementLevels, PETSC_NULL);
  CHKERRQ(ierr);
  refineIntegration = PETSC_TRUE;
  ierr =
      PetscOptionsBool("-se_refined_integration",
                       "if set element is subdivided to generate quadrature",
                       "", refineIntegration, &refineIntegration, PETSC_NULL);
  CHKERRQ(ierr);
  ierr = PetscOptionsEnd();
  CHKERRQ(ierr);
  MoFEMFunctionReturnHot(0);
}
 
// specialization
template <>
template <>
int CrackFrontSingularBase<NonlinearElasticElement::MyVolumeFE,
                           VolumeElementForcesAndSourcesCore>::
    getRuleImpl(TwoType<NonlinearElasticElement::MyVolumeFE,
                        VolumeElementForcesAndSourcesCore>,
                int order);
 
// specialization
template <>
template <>
MoFEMErrorCode CrackFrontSingularBase<NonlinearElasticElement::MyVolumeFE,
                                      VolumeElementForcesAndSourcesCore>::
    setGaussPtsImpl(TwoType<NonlinearElasticElement::MyVolumeFE,
                            VolumeElementForcesAndSourcesCore>,
                    int order);
 
typedef CrackFrontSingularBase<NonlinearElasticElement::MyVolumeFE,
                               VolumeElementForcesAndSourcesCore>
    CrackFrontElement;
 
struct OpAnalyticalSpatialTraction
    : public FaceElementForcesAndSourcesCore::UserDataOperator {
 
  const Range tRis;
  boost::ptr_vector<MethodForForceScaling> &methodsOp;
  boost::ptr_vector<NeumannForcesSurface::MethodForAnalyticalForce>
      &analyticalForceOp;
 
  CrackFrontSingularBase<VolumeElementForcesAndSourcesCoreOnSide,
                         VolumeElementForcesAndSourcesCore>
      volSideFe;
  // boost::shared_ptr<MatrixDouble> hG; ///< spatial deformation gradient
  // boost::shared_ptr<MatrixDouble> HG; ///< spatial deformation gradient
  VectorDouble sNrm; ///< Length of the normal vector
  bool setSingularCoordinates;
 
  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);
 
  VectorDouble Nf; //< Local force vector
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
};
 
struct FirendVolumeOnSide : public VolumeElementForcesAndSourcesCoreOnSide {
  using VolumeElementForcesAndSourcesCoreOnSide::
      VolumeElementForcesAndSourcesCoreOnSide;
  using UserDataOperator =
      VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator;
  using VolumeElementForcesAndSourcesCoreOnSide::
      calHierarchicalBaseFunctionsOnElement;
};
 
struct OpAnalyticalMaterialTraction
    : public FaceElementForcesAndSourcesCore::UserDataOperator {
 
  const Range tRis;
  boost::ptr_vector<MethodForForceScaling> &methodsOp;
  boost::ptr_vector<NeumannForcesSurface::MethodForAnalyticalForce>
      &analyticalForceOp;
 
  CrackFrontSingularBase<FirendVolumeOnSide, VolumeElementForcesAndSourcesCore>
      volSideFe;
  boost::shared_ptr<MatrixDouble> hG; ///< spatial deformation gradient
  boost::shared_ptr<MatrixDouble> HG; ///< spatial deformation gradient
  VectorDouble sNrm;                  ///< Length of the normal vector
  bool setSingularCoordinates;
  Range forcesOnlyOnEntitiesRow;
 
  VectorInt iNdices;
 
  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 = NULL);
 
  VectorDouble Nf; //< Local force vector
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
};
 
struct OpGetCrackFrontDataGradientAtGaussPts
    : public OpCalculateVectorFieldGradient<3, 3> {
 
  bool &singularElement;
  MatrixDouble &invSJac;
 
  OpGetCrackFrontDataGradientAtGaussPts(
      const std::string field_name, boost::shared_ptr<MatrixDouble> data_at_pts,
      bool &singular_element, MatrixDouble &inv_s_jac)
      : OpCalculateVectorFieldGradient<3, 3>(field_name, data_at_pts),
        singularElement(singular_element), invSJac(inv_s_jac) {}
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
};
 
struct OpGetCrackFrontDataAtGaussPts
    : public NonlinearElasticElement::OpGetDataAtGaussPts {
 
  bool &singularElement;
  MatrixDouble &invSJac;
 
  OpGetCrackFrontDataAtGaussPts(
      const std::string field_name,
      std::vector<VectorDouble> &values_at_gauss_pts,
      std::vector<MatrixDouble> &gradient_at_gauss_pts, bool &singular_element,
      MatrixDouble &inv_s_jac)
      : NonlinearElasticElement::OpGetDataAtGaussPts(
            field_name, values_at_gauss_pts, gradient_at_gauss_pts),
        singularElement(singular_element), invSJac(inv_s_jac) {}
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
};
 
struct OpGetCrackFrontCommonDataAtGaussPts
    : public OpGetCrackFrontDataAtGaussPts {
  OpGetCrackFrontCommonDataAtGaussPts(
      const std::string field_name,
      NonlinearElasticElement::CommonData &common_data, bool &singular_element,
      MatrixDouble &inv_s_jac)
      : OpGetCrackFrontDataAtGaussPts(field_name,
                                      common_data.dataAtGaussPts[field_name],
                                      common_data.gradAtGaussPts[field_name],
                                      singular_element, inv_s_jac) {}
};
 
struct OpPostProcDisplacements
    : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
  bool &singularElement;
  MatrixDouble &singularDisp;
  moab::Interface &postProcMesh;
  std::vector<EntityHandle> &mapGaussPts;
  boost::shared_ptr<MatrixDouble> H;
  OpPostProcDisplacements(bool &singular_element,
                          MatrixDouble &singular_ref_coords,
                          moab::Interface &post_proc_mesh,
                          std::vector<EntityHandle> &map_gauss_pts,
                          boost::shared_ptr<MatrixDouble> &H)
      : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
            "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
        singularElement(singular_element), singularDisp(singular_ref_coords),
        postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts), H(H) {}
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
};
 
struct OpTransfromSingularBaseFunctions
    : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
  bool &singularElement;
  VectorDouble &detS;
  MatrixDouble &invSJac;
  const FieldApproximationBase bAse;
  const bool applyDet;
  OpTransfromSingularBaseFunctions(
      bool &singular_element, VectorDouble &det_s, MatrixDouble &inv_s_jac,
      FieldApproximationBase base = AINSWORTH_LOBATTO_BASE,
      bool apply_det = true)
      : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(H1),
        singularElement(singular_element), detS(det_s), invSJac(inv_s_jac),
        bAse(base), applyDet(apply_det) {}
  MatrixDouble diffNinvJac;
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
};
 
/**
* \brief Calculate explicit derivative of energy
 
\f[
\left(
\frac{\partial \Psi}{\partial \mathbf{X}}
\right)_\textrm{expl.}=
\left.\left(
\frac{\partial \Psi}{\partial \mathbf{X}}
\right)\right|_{\mathbf{u}=const.,\mathbf{F}=const.}
\f]
]
 
*/
struct OpRhsBoneExplicitDerivariveWithHooke
    : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
  HookeElement::DataAtIntegrationPts &commonData;
  boost::shared_ptr<VectorDouble> rhoAtGaussPts;
  boost::shared_ptr<MatrixDouble> rhoGradAtGaussPts;
  Range tetsInBlock;
  const double rHo0;
  const double boneN;
  Range forcesOnlyOnEntitiesRow;
  ublas::vector<int> iNdices;
 
  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 = NULL);
  VectorDouble nF;
  MoFEMErrorCode doWork(int row_side, EntityType row_type,
                        DataForcesAndSourcesCore::EntData &row_data);
};
/**
* \brief Calculate explicit derivative of energy
 
\f[
\left(
\frac{\partial \Psi}{\partial \mathbf{X}}
\right)_\textrm{expl.}=
\left.\left(
\frac{\partial \Psi}{\partial \mathbf{X}}
\right)\right|_{\mathbf{u}=const.,\mathbf{F}=const.}
\f]
]
 
*/
struct OpLhsBoneExplicitDerivariveWithHooke_dX
    : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
  HookeElement::DataAtIntegrationPts &commonData;
  boost::shared_ptr<VectorDouble> rhoAtGaussPts;
  boost::shared_ptr<MatrixDouble> diffRhoAtGaussPts;
  boost::shared_ptr<MatrixDouble> diffDiffRhoAtGaussPts;
  MatrixDouble &singularDispl;
  Range tetsInBlock;
  const double rHo0;
  const double boneN;
  Range forcesOnlyOnEntitiesRow;
  ublas::vector<int> iNdices;
  MatrixDouble nA;
  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 = NULL);
  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                        EntityType col_type,
                        DataForcesAndSourcesCore::EntData &row_data,
                        DataForcesAndSourcesCore::EntData &col_data);
};
 
struct OpLhsBoneExplicitDerivariveWithHooke_dx
    : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
  HookeElement::DataAtIntegrationPts &commonData;
  boost::shared_ptr<VectorDouble> rhoAtGaussPts;
  boost::shared_ptr<MatrixDouble> diffRhoAtGaussPts;
  boost::shared_ptr<MatrixDouble> diffDiffRhoAtGaussPts;
  MatrixDouble &singularDispl;
  Range tetsInBlock;
  const double rHo0;
  const double boneN;
  Range forcesOnlyOnEntitiesRow;
  ublas::vector<int> iNdices;
  MatrixDouble nA;
  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 = NULL);
  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                        EntityType col_type,
                        DataForcesAndSourcesCore::EntData &row_data,
                        DataForcesAndSourcesCore::EntData &col_data);
};
 
struct OpAleLhsWithDensitySingularElement_dX_dX
    : public VolumeElementForcesAndSourcesCore::UserDataOperator {
 
  boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
  boost::shared_ptr<MatrixDouble> rhoGradAtGaussPtsPtr;
  boost::shared_ptr<MatrixDouble> singularDisplacement;
 
  const double rhoN;
  const double rHo0;
 
  // Finite element stiffness sub-matrix K_ij
  MatrixDouble K;
  MatrixDouble transK;
  VectorDouble nF;
 
  boost::shared_ptr<HookeElement::DataAtIntegrationPts> dataAtPts;
 
  VectorInt rowIndices;
  VectorInt colIndices;
 
  int nbRows;           ///< number of dofs on rows
  int nbCols;           ///< number if dof on column
  int nbIntegrationPts; ///< number of integration points
  bool isDiag;          ///< true if this block is on diagonal
 
  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);
 
protected:
  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                        EntityType col_type,
                        DataForcesAndSourcesCore::EntData &row_data,
                        DataForcesAndSourcesCore::EntData &col_data);
 
  MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data,
                           DataForcesAndSourcesCore::EntData &col_data);
 
  MoFEMErrorCode aSsemble(DataForcesAndSourcesCore::EntData &row_data,
                          DataForcesAndSourcesCore::EntData &col_data);
};
 
struct OpAleLhsWithDensitySingularElement_dx_dX
    : public VolumeElementForcesAndSourcesCore::UserDataOperator {
 
  boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
  boost::shared_ptr<MatrixDouble> rhoGradAtGaussPtsPtr;
  boost::shared_ptr<MatrixDouble> singularDisplacement;
 
  const double rhoN;
  const double rHo0;
 
  // Finite element stiffness sub-matrix K_ij
  MatrixDouble K;
  MatrixDouble transK;
  VectorDouble nF;
 
  boost::shared_ptr<HookeElement::DataAtIntegrationPts> dataAtPts;
 
  VectorInt rowIndices;
  VectorInt colIndices;
 
  int nbRows;           ///< number of dofs on rows
  int nbCols;           ///< number if dof on column
  int nbIntegrationPts; ///< number of integration points
  bool isDiag;          ///< true if this block is on diagonal
 
  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);
 
protected:
  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                        EntityType col_type,
                        DataForcesAndSourcesCore::EntData &row_data,
                        DataForcesAndSourcesCore::EntData &col_data);
 
  MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data,
                           DataForcesAndSourcesCore::EntData &col_data);
 
  MoFEMErrorCode aSsemble(DataForcesAndSourcesCore::EntData &row_data,
                          DataForcesAndSourcesCore::EntData &col_data);
};
 
/**
 * \brief Op to generate artificial density field
 *
 */
struct OpGetDensityFieldForTesting : public HookeElement::VolUserDataOperator {
 
  boost::shared_ptr<MatrixDouble> matCoordsPtr;
  boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
  boost::shared_ptr<MatrixDouble> rhoGradAtGaussPtsPtr;
  boost::shared_ptr<MatrixDouble> rhoGradGradAtGaussPtsPtr;
  boost::shared_ptr<CrackFrontElement> feSingularPtr;
  MatrixDouble &singularDisp;
  boost::shared_ptr<MatrixDouble> matGradPosAtPtsPtr;
  OpGetDensityFieldForTesting(
      const std::string row_field,
      boost::shared_ptr<MatrixDouble> mat_coords_ptr,
      boost::shared_ptr<VectorDouble> density_at_pts,
      boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts_ptr,
      boost::shared_ptr<MatrixDouble> rho_grad_grad_at_gauss_pts_ptr,
      boost::shared_ptr<CrackFrontElement> fe_singular_ptr,
      MatrixDouble &singular_disp,
      boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr)
      : HookeElement::VolUserDataOperator(row_field, OPROW),
        matCoordsPtr(mat_coords_ptr), rhoAtGaussPtsPtr(density_at_pts),
        rhoGradAtGaussPtsPtr(rho_grad_at_gauss_pts_ptr),
        rhoGradGradAtGaussPtsPtr(rho_grad_grad_at_gauss_pts_ptr),
        feSingularPtr(fe_singular_ptr), singularDisp(singular_disp),
        matGradPosAtPtsPtr(mat_grad_pos_at_pts_ptr) {}
 
  MoFEMErrorCode doWork(int row_side, EntityType row_type,
                        HookeElement::EntData &row_data);
};
 
/**
 * \brief Mark crack surfaces on skin
 *
 */
struct OpSetTagRangeOnSkin
    : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
 
  moab::Interface &postProcMesh;
  std::vector<EntityHandle> &mapGaussPts;
  const string tagName;
  Range rangeToTag;
  int tagValue;
 
  OpSetTagRangeOnSkin(moab::Interface &post_proc_mesh,
                      vector<EntityHandle> &map_gauss_pts, Range &range_to_tag,
                      const string &tag_name, int tag_value = 1)
      : FaceElementForcesAndSourcesCore::UserDataOperator(
            "MESH_NODE_POSITIONS", UserDataOperator::OPCOL),
        postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
        rangeToTag(range_to_tag), tagName(tag_name), tagValue(tag_value) {}
 
  MoFEMErrorCode doWork(int row_side, EntityType row_type,
                        DataForcesAndSourcesCore::EntData &row_data);
};
 
} // namespace FractureMechanics
 
#endif // __CRACK_FRONT_ELEMENT_HPP__