GriffithForceElement.hpp

Go to the documentation of this file.

/* 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 __GRIFFITH_FORCE_ELEMENT_HPP__
#define __GRIFFITH_FORCE_ELEMENT_HPP__
 
#ifndef WITH_ADOL_C
#error "MoFEM need to be compiled with ADOL-C"
#endif
 
namespace FractureMechanics {
 
static double calMax(double a, double b, double r) {
  return (a + b + (1 / r) * pow(fabs(a - b), r)) / 2;
}
 
static double diffCalMax_a(double a, double b, double r) {
  double sgn = ((a - b) == 0) ? 0 : (((a - b) < 0) ? -1 : 1);
  return (1 + pow(fabs(a - b), r - 1) * sgn * (+1)) / 2;
}
 
/** \brief Implementation of Griffith element
 */
struct GriffithForceElement {
 
  MoFEM::Interface &mField;
 
  struct MyTriangleFE : public FaceElementForcesAndSourcesCore {
 
    Mat B;
    Vec F;
 
    MyTriangleFE(MoFEM::Interface &m_field)
        : FaceElementForcesAndSourcesCore(m_field), B(PETSC_NULL),
          F(PETSC_NULL) {}
 
    int getRule(int order) { return 1; };
  };
 
  boost::shared_ptr<MyTriangleFE> feRhsPtr;
  boost::shared_ptr<MyTriangleFE> feLhsPtr;
 
  MyTriangleFE &feRhs;
  MyTriangleFE &getLoopFeRhs() { return feRhs; }
  MyTriangleFE &feLhs;
  MyTriangleFE &getLoopFeLhs() { return feLhs; }
 
  GriffithForceElement(MoFEM::Interface &m_field)
      : mField(m_field), feRhsPtr(new MyTriangleFE(m_field)),
        feLhsPtr(new MyTriangleFE(m_field)), feRhs(*feRhsPtr),
        feLhs(*feLhsPtr) {}
 
  struct CommonData {
    VectorDouble griffithForce;        ///< Griffith force at nodes
    MatrixDouble tangentGriffithForce; ///< Tangent matrix
    vector<double *>
        tangentGriffithForceRowPtr; ///< Pointers to rows of tangent matrix
    VectorDouble densityRho;        ///< gc * rho^beta
    MatrixDouble diffRho;           ///< for lhs with heterogeneous gc
  };
  CommonData commonData; ///< Common data sgared between operators
 
  struct BlockData {
    double gc;        ///< Griffith energy
    double E;         ///< Young modulus
    double r;         ///< regularity parameter
    Range frontNodes; ///< Nodes on crack front
    BlockData() : r(1) {}
  };
  map<int, BlockData> blockData; ///< Shared data between finite elements
 
  static MoFEMErrorCode getElementOptions(BlockData &block_data) {
    MoFEMFunctionBegin;
 
    MOFEM_LOG_CHANNEL("WORLD");
    MOFEM_LOG_FUNCTION();
    MOFEM_LOG_TAG("WORLD", "Griffith");
 
    ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "",
                             "Get Griffith element options", "none");
    CHKERRQ(ierr);
    CHKERR PetscOptionsScalar("-griffith_gc", "release energy", "",
                              block_data.gc, &block_data.gc, PETSC_NULL);
    CHKERR PetscOptionsScalar("-griffith_r", "release energy", "", block_data.r,
                              &block_data.r, PETSC_NULL);
    CHKERR PetscOptionsScalar("-griffith_E", "stiffness", "", block_data.E,
                              &block_data.E, PETSC_NULL);
    ierr = PetscOptionsEnd();
    CHKERRQ(ierr);
 
    MOFEM_LOG_C("WORLD", Sev::inform, "### Input parameter: -griffith_E %6.4e",
                block_data.E);
    MOFEM_LOG_C("WORLD", Sev::inform, "### Input parameter: -griffith_r %6.4e",
                block_data.r);
 
    MoFEMFunctionReturn(0);
  }
 
  /**
 
      \f[
      A = \| \mathbf{N} \| =
        \left\|
          \textrm{Spin}\left[
            \frac{\partial\mathbf{N}\mathbf{X}}{\partial\xi}
          \right]
          \frac{\partial\mathbf{N}\mathbf{X}}{\partial\eta}
          \right\|
      \f]
 
      */
  template <typename TYPE> struct AuxFunctions {
 
    ublas::vector<TYPE> referenceCoords;
    ublas::vector<TYPE> currentCoords;
    ublas::vector<TYPE> griffithForce;
 
    FTensor::Tensor1<TYPE, 3> t_Normal;
    FTensor::Tensor1<TYPE, 3> t_Tangent1_ksi;
    FTensor::Tensor1<TYPE, 3> t_Tangent2_eta;
    TYPE currentArea;
 
    FTensor::Index<'i', 3> i;
    FTensor::Index<'j', 3> j;
    FTensor::Index<'k', 3> k;
    FTensor::Number<0> N0;
    FTensor::Number<1> N1;
 
    boost::function<FTensor::Tensor1<TYPE *, 3>(ublas::vector<TYPE> &v)>
        get_ftensor_from_vec;
 
    boost::function<FTensor::Tensor1<FTensor::PackPtr<const double *, 2>, 2>(
        const MatrixDouble &m)>
        get_ftensor_from_mat_2d;
 
    AuxFunctions() {
 
      get_ftensor_from_vec = [](ublas::vector<TYPE> &v) {
        return FTensor::Tensor1<TYPE *, 3>(&v(0), &v(1), &v(2), 3);
      };
 
      get_ftensor_from_mat_2d = [](const MatrixDouble &m) {
        return FTensor::Tensor1<FTensor::PackPtr<const double *, 2>, 2>(
            &m(0, 0), &m(0, 1));
      };
    }
 
    MoFEMErrorCode calculateAreaAndNormal(const MatrixDouble &diff_n) {
      MoFEMFunctionBeginHot;
 
      auto t_coords = get_ftensor_from_vec(currentCoords);
      auto t_diff = get_ftensor_from_mat_2d(diff_n);
      t_Tangent1_ksi(i) = (TYPE)0;
      t_Tangent2_eta(i) = (TYPE)0;
      for (int nn = 0; nn != 3; ++nn) {
        t_Tangent1_ksi(i) += t_coords(i) * t_diff(N0);
        t_Tangent2_eta(i) += t_coords(i) * t_diff(N1);
        ++t_diff;
        ++t_coords;
      }
 
      t_Normal(j) =
          FTensor::levi_civita(i, j, k) * t_Tangent1_ksi(k) * t_Tangent2_eta(i);
      currentArea = 0.5 * sqrt(t_Normal(i) * t_Normal(i));
 
      MoFEMFunctionReturnHot(0);
    }
 
    MoFEMErrorCode calculateMatA(const MatrixDouble &diff_n) {
      MoFEMFunctionBeginHot;
 
      griffithForce.resize(9, false);
      std::fill(griffithForce.begin(), griffithForce.end(), 0);
      auto t_griffith = get_ftensor_from_vec(griffithForce);
      auto t_diff = get_ftensor_from_mat_2d(diff_n);
 
      for (int dd = 0; dd != 3; dd++) {
 
        t_griffith(i) += FTensor::levi_civita(i, j, k) * t_Tangent1_ksi(k) *
                         t_diff(N1) * t_Normal(j) * 0.25 / currentArea;
        t_griffith(i) -= FTensor::levi_civita(i, j, k) * t_Tangent2_eta(k) *
                         t_diff(N0) * t_Normal(j) * 0.25 / currentArea;
        ++t_diff;
        ++t_griffith;
      }
 
      MoFEMFunctionReturnHot(0);
    }
 
    MoFEMErrorCode calculateGriffithForce(const double gc, const double beta,
                                          const MatrixDouble &diff_n) {
      MoFEMFunctionBeginHot;
 
      CHKERR calculateAreaAndNormal(diff_n);
      CHKERR calculateMatA(diff_n);
 
      for (int dd = 0; dd != 9; dd++)
        griffithForce[dd] *= gc * beta;
 
      MoFEMFunctionReturnHot(0);
    }
  };
 
  struct AuxOp {
 
    int tAg;
    BlockData &blockData;
    CommonData &commonData;
 
    AuxOp(int tag, BlockData &block_data, CommonData &common_data)
        : tAg(tag), blockData(block_data), commonData(common_data){};
 
    ublas::vector<int> rowIndices;
    ublas::vector<int> rowIndicesLocal;
 
    VectorInt3 rowLambdaIndices;
    VectorInt3 rowLambdaIndicesLocal;
    VectorDouble3 lambdaAtNodes;
    VectorDouble activeVariables;
 
    MoFEMErrorCode setIndices(DataForcesAndSourcesCore::EntData &data) {
      MoFEMFunctionBegin;
      const auto nb_dofs = data.getIndices().size();
      rowIndices.resize(nb_dofs, false);
      noalias(rowIndices) = data.getIndices();
      rowIndicesLocal.resize(nb_dofs, false);
      noalias(rowIndicesLocal) = data.getLocalIndices();
      auto dit = data.getFieldDofs().begin();
      auto hi_dit = data.getFieldDofs().end();
      for (int ii = 0; dit != hi_dit; ++dit, ++ii) {
        if (auto dof = (*dit)) {
          if (blockData.frontNodes.find(dof->getEnt()) ==
              blockData.frontNodes.end()) {
 
            const auto side = dof->getSideNumber();
            const auto idx = dof->getEntDofIdx();
            const auto rank = dof->getNbDofsOnEnt();
            if (ii != side * rank + idx)
              SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong index");
 
            rowIndices[ii] = -1;
            rowIndicesLocal[ii] = -1;
          }
        }
      }
      MoFEMFunctionReturn(0);
    }
 
    MoFEMErrorCode
    setVariables(FaceElementForcesAndSourcesCore::UserDataOperator *fe_ptr,
                 DataForcesAndSourcesCore::EntData &data) {
      MoFEMFunctionBeginHot;
      int nb_dofs = data.getIndices().size();
      activeVariables.resize(18, false);
      for (int dd = 0; dd != nb_dofs; dd++) {
        activeVariables[dd] = fe_ptr->getCoords()[dd];
      }
      for (int dd = 0; dd != nb_dofs; dd++) {
        activeVariables[9 + dd] = data.getFieldData()[dd];
      }
      MoFEMFunctionReturnHot(0);
    }
 
    MoFEMErrorCode
    setLambdaNodes(FaceElementForcesAndSourcesCore::UserDataOperator *fe_ptr,
                   const std::string &lambda_field_name) {
      MoFEMFunctionBeginHot;
      lambdaAtNodes.resize(3, false);
      lambdaAtNodes.clear();
      const auto bit_field_number =
          fe_ptr->getFEMethod()->getFieldBitNumber(lambda_field_name);
 
      auto data_dofs = fe_ptr->getFEMethod()->getDataDofsPtr();
      for (auto dit = data_dofs->get<Unique_mi_tag>().lower_bound(
               FieldEntity::getLoBitNumberUId(bit_field_number));
           dit != data_dofs->get<Unique_mi_tag>().upper_bound(
                      FieldEntity::getHiBitNumberUId(bit_field_number));
           ++dit) {
        if (dit->get()->getEntType() != MBVERTEX) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "DOFs only on vertices");
        }
        int side = dit->get()->getSideNumber();
        lambdaAtNodes[side] = dit->get()->getFieldData();
      }
      MoFEMFunctionReturnHot(0);
    }
 
    MoFEMErrorCode
    setLambdaIndices(FaceElementForcesAndSourcesCore::UserDataOperator *fe_ptr,
                     const std::string &lambda_field_name) {
      MoFEMFunctionBeginHot;
      rowLambdaIndices.resize(3, false);
      rowLambdaIndicesLocal.resize(3, false);
      std::fill(rowLambdaIndices.begin(), rowLambdaIndices.end(), -1);
      std::fill(rowLambdaIndicesLocal.begin(), rowLambdaIndicesLocal.end(), -1);
      const auto bit_field_number =
          fe_ptr->getFEMethod()->getFieldBitNumber(lambda_field_name);
 
      auto data_dofs = fe_ptr->getFEMethod()->getRowDofsPtr();
      for (auto dit = data_dofs->get<Unique_mi_tag>().lower_bound(
               FieldEntity::getLoBitNumberUId(bit_field_number));
           dit != data_dofs->get<Unique_mi_tag>().upper_bound(
                      FieldEntity::getHiBitNumberUId(bit_field_number));
           ++dit) {
 
        if (dit->get()->getEntType() != MBVERTEX)
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "DOFs only on vertices");
 
        int side = dit->get()->getSideNumber();
        rowLambdaIndices[side] = dit->get()->getPetscGlobalDofIdx();
        rowLambdaIndicesLocal[side] = dit->get()->getPetscLocalDofIdx();
      }
      MoFEMFunctionReturnHot(0);
    }
  };
 
  struct OpRhs : public FaceElementForcesAndSourcesCore::UserDataOperator,
                 AuxOp {
    OpRhs(int tag, BlockData &block_data, CommonData &common_data)
        : FaceElementForcesAndSourcesCore::UserDataOperator(
              "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
          AuxOp(tag, block_data, common_data) {}
    MoFEMErrorCode doWork(int side, EntityType type,
                          DataForcesAndSourcesCore::EntData &data) {
      MoFEMFunctionBegin;
 
      if (type != MBVERTEX)
        MoFEMFunctionReturnHot(0);
      CHKERR setIndices(data);
      CHKERR setVariables(this, data);
      int nb_dofs = data.getIndices().size();
 
      CHKERR VecSetValues(
          getFEMethod()->snes_f, nb_dofs, &*rowIndices.data().begin(),
          &*commonData.griffithForce.data().begin(), ADD_VALUES);
      MoFEMFunctionReturn(0);
    }
  };
 
  /**
   * @brief calculate griffith force vector
   *
   *  \f[
       \mathbf{f}_{\mathrm{m}}^{\mathrm{h}}=\frac{1}{2}\left(\tilde{\mathbf{A}}_{\Gamma}^{\mathrm{h}}\right)^{\mathrm{T}}
  \mathbf{g}_{c}(\rho(\mathbf{X}))
   *  ]\f
   *
   * The matrix $\mathbf{A}_{\Gamma}^h$ comprising of direction vectors along
  the crack front that are normal to the crack front and tangent to the crack
  surface is defined as follows:
  \f[
    A^h_{\Gamma} =
          \| \mathbf{N}(\tilde{\mathbf{X}}) \| =
  \left\|
  \epsilon_{ijk}
  \frac{\partial \Phi^\alpha_p}{\partial \xi_i}
  \frac{\partial \Phi^\beta_r}{\partial \xi_j}
  \tilde{X}^\alpha_p
  \tilde{X}^\beta_r
  \right\|
    \f]
   *
   *
   */
  struct OpCalculateGriffithForce
      : public FaceElementForcesAndSourcesCore::UserDataOperator,
        AuxOp {
    OpCalculateGriffithForce(int tag, BlockData &block_data,
                             CommonData &common_data)
        : FaceElementForcesAndSourcesCore::UserDataOperator(
              "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
          AuxOp(tag, block_data, common_data) {}
    AuxFunctions<double> auxFun;
    MoFEMErrorCode doWork(int side, EntityType type,
                          DataForcesAndSourcesCore::EntData &data) {
      MoFEMFunctionBegin;
 
      if (type != MBVERTEX)
        MoFEMFunctionReturnHot(0);
      CHKERR setIndices(data);
      CHKERR setVariables(this, data);
      int nb_dofs = data.getIndices().size();
      int nb_gauss_pts = data.getN().size1();
 
      auxFun.referenceCoords.resize(nb_dofs, false);
      auxFun.currentCoords.resize(nb_dofs, false);
 
      for (int dd = 0; dd != nb_dofs; dd++) {
        auxFun.referenceCoords[dd] = getCoords()[dd];
        auxFun.currentCoords[dd] = data.getFieldData()[dd];
      }
 
      for (int gg = 0; gg != nb_gauss_pts; gg++) {
        double val = getGaussPts()(2, gg) * 0.5;
        CHKERR auxFun.calculateGriffithForce(blockData.gc, val,
                                             data.getDiffN(gg));
      }
      commonData.griffithForce.resize(nb_dofs, false);
      std::fill(commonData.griffithForce.begin(),
                commonData.griffithForce.end(), 0);
      std::copy(auxFun.griffithForce.begin(), auxFun.griffithForce.end(),
                commonData.griffithForce.begin());
 
      MoFEMFunctionReturn(0);
    }
  };
 
  /**
   * \brief Calculate density at the crack fron nodes and multiplity gc
   *
   *  For heterogeneous case, gc is multiplied by a coefficient depending on
   spatial distribution of density as follows:
   *
      \f[
          g_c(\mathbf X) = g_c \cdot \rho(\mathbf X)^{\beta^{gc}}
      \f]
 
  */
  struct OpCalculateRhoAtCrackFrontUsingPrecalulatedCoeffs
      : public FaceElementForcesAndSourcesCore::UserDataOperator,
        AuxOp {
    const double betaGC;
    std::string rhoTagName;
    boost::shared_ptr<MWLSApprox> mwlsApprox;
    MoFEM::Interface &mField;
 
    OpCalculateRhoAtCrackFrontUsingPrecalulatedCoeffs(
        string rho_tag_name, const double beta_gc,
        boost::shared_ptr<MWLSApprox> mwls_approx, MoFEM::Interface &m_field,
        int tag, BlockData &block_data, CommonData &common_data)
        : FaceElementForcesAndSourcesCore::UserDataOperator(
              "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
          rhoTagName(rho_tag_name), betaGC(beta_gc), mwlsApprox(mwls_approx),
          mField(m_field), AuxOp(tag, block_data, common_data) {}
 
    MoFEMErrorCode doWork(int side, EntityType type,
                          DataForcesAndSourcesCore::EntData &row_data) {
      MoFEMFunctionBegin;
      FTensor::Index<'i', 3> i;
      FTensor::Index<'j', 3> j;
 
      if (type != MBVERTEX)
        MoFEMFunctionReturnHot(0);
      const int row_nb_dofs = row_data.getIndices().size();
      if (!row_nb_dofs)
        MoFEMFunctionReturnHot(0);
 
      CHKERR setIndices(row_data);
      CHKERR setVariables(this, row_data);
 
      auto &inv_AB_map = mwlsApprox->invABMap.at(getFEEntityHandle());
      auto &influence_nodes_map =
          mwlsApprox->influenceNodesMap.at(getFEEntityHandle());
      auto &dm_nodes_map = mwlsApprox->dmNodesMap[this->getFEEntityHandle()];
 
      auto get_ftensor_from_vec = [](auto &v) {
        return FTensor::Tensor1<double *, 3>(&v(0), &v(1), &v(2), 3);
      };
 
      Tag th;
      CHKERR mwlsApprox->mwlsMoab.tag_get_handle(rhoTagName.c_str(), th);
 
      int nb_dofs = row_data.getFieldData().size();
 
      commonData.diffRho.resize(3, nb_dofs / 3, false);
      std::fill(commonData.diffRho.data().begin(),
                commonData.diffRho.data().end(), 0);
 
      auto t_coords = get_ftensor_from_vec(row_data.getFieldData());
 
      commonData.densityRho.resize(nb_dofs / 3, false);
      std::fill(commonData.densityRho.begin(), commonData.densityRho.end(), 0);
      auto rho = getFTensor0FromVec(commonData.densityRho);
      VectorDouble gc_coeff(nb_dofs, false);
      auto t_coeff = get_ftensor_from_vec(gc_coeff);
      auto t_approx_base_point =
          getFTensor1FromMat<3>(mwlsApprox->approxBasePoint);
      const int nb_gauss_pts = row_data.getN().size1();
 
      for (int gg = 0; gg != nb_gauss_pts; ++gg) {
        for (int dd = 0; dd != nb_dofs / 3; ++dd) {
 
          // go only through nodes at the crack front
          if (rowIndices[dd * 3] != -1) {
            FTensor::Tensor1<double, 3> t_pos;
            t_pos(i) = t_coords(i);
 
            for (int dd : {0, 1, 2})
              mwlsApprox->approxPointCoords[dd] = t_approx_base_point(dd);
 
            mwlsApprox->getInfluenceNodes() = influence_nodes_map[gg];
            mwlsApprox->getInvAB() = inv_AB_map[gg];
            mwlsApprox->getShortenDm() = dm_nodes_map[gg];
 
            // approximate at MESH_NODE_POSITIONS
            CHKERR mwlsApprox->evaluateApproxFun(&t_pos(0));
            CHKERR mwlsApprox->evaluateFirstDiffApproxFun(&t_pos(0), false);
            CHKERR mwlsApprox->getTagData(th);
            rho = mwlsApprox->getDataApprox()[0];
            const auto &diff_vals = mwlsApprox->getDiffDataApprox();
            for (int ii = 0; ii != 3; ++ii) {
              commonData.diffRho(ii, dd) = diff_vals(ii, 0);
            }
            t_coeff(i) = pow(rho, betaGC);
          }
 
          ++rho;
          ++t_coords;
          ++t_coeff;
        }
      }
      // multiply griffith force with coefficients
      commonData.griffithForce =
          element_prod(commonData.griffithForce, gc_coeff);
      MoFEMFunctionReturn(0);
    }
  };
 
  struct SaveGriffithForceOnTags {
 
    MoFEM::Interface &mField;
 
    Tag tH;
    Range &tRis;
    Tag *thCoordsPtr;
 
    SaveGriffithForceOnTags(MoFEM::Interface &m_field, Tag th, Range &tris,
                            Tag *th_coords = nullptr)
        : mField(m_field), tH(th), tRis(tris), thCoordsPtr(th_coords) {}
 
    AuxFunctions<double> auxFun;
 
    MoFEMErrorCode operator()() {
      MoFEMFunctionBegin;
 
      Range verts;
      CHKERR mField.get_moab().get_connectivity(tRis, verts, true);
      {
        std::vector<double> data(verts.size() * 3, 0);
        CHKERR mField.get_moab().tag_set_data(tH, verts, &*data.begin());
      }
      MatrixDouble diff_n(3, 2);
      CHKERR ShapeDiffMBTRI(&*diff_n.data().begin());
 
      for (Range::iterator tit = tRis.begin(); tit != tRis.end(); ++tit) {
        int num_nodes;
        const EntityHandle *conn;
        CHKERR mField.get_moab().get_connectivity(*tit, conn, num_nodes, true);
        VectorDouble9 coords(9);
 
        int nb_dofs = num_nodes * 3;
        auxFun.referenceCoords.resize(9, false);
        auxFun.currentCoords.resize(9, false);
 
        if (thCoordsPtr) {
          CHKERR mField.get_moab().tag_get_data(*thCoordsPtr, conn, num_nodes,
                                                &*coords.begin());
        } else {
          CHKERR mField.get_moab().get_coords(conn, num_nodes,
                                              &*coords.begin());
        }
        for (int dd = 0; dd != nb_dofs; dd++) {
          auxFun.currentCoords[dd] = coords[dd];
        }
        if (thCoordsPtr)
          CHKERR mField.get_moab().get_coords(conn, num_nodes,
                                              &*coords.begin());
        for (int dd = 0; dd != nb_dofs; dd++) {
          auxFun.referenceCoords[dd] = coords[dd];
        }
 
        CHKERR auxFun.calculateGriffithForce(1, 0.5, diff_n);
 
        double *tag_vals[3];
        CHKERR mField.get_moab().tag_get_by_ptr(tH, conn, 3,
                                                (const void **)tag_vals);
        for (int nn = 0; nn != 3; ++nn) {
          for (int dd = 0; dd != 3; ++dd) {
            tag_vals[nn][dd] += auxFun.griffithForce[3 * nn + dd];
          }
        }
      }
 
      MoFEMFunctionReturn(0);
    }
  };
 
  struct OpLhs : public FaceElementForcesAndSourcesCore::UserDataOperator,
                 AuxOp {
 
    OpLhs(int tag, BlockData &block_data, CommonData &common_data,
          const double beta_gc = 0)
        : FaceElementForcesAndSourcesCore::UserDataOperator(
              "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
              UserDataOperator::OPROWCOL),
          betaGC(beta_gc), AuxOp(tag, block_data, common_data) {}
 
    MatrixDouble mat;
    const double betaGC;
 
    MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                          EntityType col_type,
                          DataForcesAndSourcesCore::EntData &row_data,
                          DataForcesAndSourcesCore::EntData &col_data) {
      MoFEMFunctionBegin;
 
      if (row_type != MBVERTEX) {
        MoFEMFunctionReturnHot(0);
      }
 
      FTensor::Index<'i', 3> i;
      FTensor::Index<'j', 3> j;
      FTensor::Index<'l', 3> l;
      FTensor::Index<'m', 3> m;
      FTensor::Index<'k', 3> k;
      FTensor::Number<0> N0;
      FTensor::Number<1> N1;
 
      AuxFunctions<double> auxFun;
      auxFun.currentCoords.resize(9, false);
 
      int nb_dofs = row_data.getFieldData().size();
      for (int dd = 0; dd != nb_dofs; dd++)
        auxFun.currentCoords[dd] = row_data.getFieldData()[dd];
 
      CHKERR setIndices(row_data);
      CHKERR setVariables(this, row_data);
      int nb_gauss_pts = row_data.getN().size1();
      int nb_base_row = row_data.getFieldData().size() / 3;
      int nb_base_col = col_data.getFieldData().size() / 3;
      int row_nb_dofs = row_data.getIndices().size();
 
      mat.resize(9, 9, false);
      mat.clear();
      auto calculate_derivative = [&](FTensor::Tensor1<double *, 2> &t_diff) {
        FTensor::Tensor2<double, 3, 3> t_d_n;
        t_d_n(i, j) = FTensor::levi_civita(i, j, k) * auxFun.t_Tangent1_ksi(k) *
                      t_diff(N1);
        t_d_n(i, j) -= FTensor::levi_civita(i, j, k) *
                       auxFun.t_Tangent2_eta(k) * t_diff(N0);
 
        return t_d_n;
      };
 
      CHKERR auxFun.calculateAreaAndNormal(col_data.getDiffN());
      auto t_normal = auxFun.t_Normal;
      const double coefficient_1 = 0.5 * pow((t_normal(i) * t_normal(i)), -0.5);
      const double coefficient_2 = 0.5 * pow((t_normal(i) * t_normal(i)), -1.5);
 
      // for homogeneous case, this is a bit ugly
      auto &rho_v = commonData.densityRho;
      if (rho_v.empty() || rho_v.size() != nb_base_row) {
        rho_v.resize(nb_base_row, false);
      }
      auto rho = getFTensor0FromVec(rho_v);
      for (int gg = 0; gg != nb_gauss_pts; gg++) {
 
        double val = getGaussPts()(2, gg) * 0.5;
        auto t_row_diff = row_data.getFTensor1DiffN<2>(gg, 0);
        for (int rrr = 0; rrr != nb_base_row; rrr++) {
 
          FTensor::Tensor2<double *, 3, 3> t_mat(
              &mat(3 * rrr + 0, 0), &mat(3 * rrr + 0, 1), &mat(3 * rrr + 0, 2),
              &mat(3 * rrr + 1, 0), &mat(3 * rrr + 1, 1), &mat(3 * rrr + 1, 2),
              &mat(3 * rrr + 2, 0), &mat(3 * rrr + 2, 1), &mat(3 * rrr + 2, 2),
              3);
 
          auto t_tan_row = calculate_derivative(t_row_diff);
          auto t_col_diff = col_data.getFTensor1DiffN<2>(gg, 0);
          const double coef = blockData.gc * val * pow(rho, betaGC);
          for (int ccc = 0; ccc != nb_base_col; ccc++) {
 
            auto t_tan_col = calculate_derivative(t_col_diff);
 
            t_mat(i, j) -= coefficient_1 *
                           (t_tan_row(i, l) * t_tan_col(l, j) +
                            FTensor::levi_civita(i, j, k) * t_col_diff(N0) *
                                t_row_diff(N1) * t_normal(k) -
                            FTensor::levi_civita(i, j, k) * t_col_diff(N1) *
                                t_row_diff(N0) * t_normal(k));
 
            t_mat(i, j) += coefficient_2 * (t_tan_row(i, k) * t_normal(k)) *
                           (t_tan_col(l, j) * t_normal(l));
 
            t_mat(i, j) *= coef;
            ++t_col_diff;
            ++t_mat;
          }
          ++t_row_diff;
          ++rho;
        }
      }
 
      int col_nb_dofs = col_data.getIndices().size();
 
      CHKERR MatSetValues(getFEMethod()->snes_B, row_nb_dofs,
                          &*rowIndices.data().begin(), col_nb_dofs,
                          &*col_data.getIndices().data().begin(),
                          &*mat.data().begin(), ADD_VALUES);
 
      MoFEMFunctionReturn(0);
    }
  };
 
  struct OpHeterogeneousGcLhs
      : public FaceElementForcesAndSourcesCore::UserDataOperator,
        AuxOp {
 
    OpHeterogeneousGcLhs(int tag, BlockData &block_data,
                         CommonData &common_data, const double beta_gc)
        : FaceElementForcesAndSourcesCore::UserDataOperator(
              "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
              UserDataOperator::OPROWCOL),
          betaGC(beta_gc), AuxOp(tag, block_data, common_data) {}
 
    MatrixDouble mat;
    const double betaGC;
 
    MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                          EntityType col_type,
                          DataForcesAndSourcesCore::EntData &row_data,
                          DataForcesAndSourcesCore::EntData &col_data) {
      MoFEMFunctionBegin;
 
      if (row_type != MBVERTEX) {
        MoFEMFunctionReturnHot(0);
      }
 
      FTensor::Index<'i', 3> i;
      FTensor::Index<'j', 3> j;
 
      AuxFunctions<double> auxFun;
      auxFun.currentCoords.resize(9, false);
 
      int nb_dofs = row_data.getFieldData().size();
      for (int dd = 0; dd != nb_dofs; dd++)
        auxFun.currentCoords[dd] = row_data.getFieldData()[dd];
 
      CHKERR setIndices(row_data);
      CHKERR setVariables(this, row_data);
      int nb_gauss_pts = row_data.getN().size1();
      int nb_base_row = row_data.getFieldData().size() / 3;
      int nb_base_col = col_data.getFieldData().size() / 3;
 
      int row_nb_dofs = row_data.getIndices().size();
      mat.resize(9, 9, false);
      mat.clear();
 
      const double &gc = blockData.gc;
      auto rho = getFTensor0FromVec(commonData.densityRho);
      auto get_ftensor_from_vec = [](auto &v) {
        return FTensor::Tensor1<double *, 3>(&v(0), &v(1), &v(2), 3);
      };
 
      auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
        return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(
            &m(3 * r + 0, 3 * c + 0), &m(3 * r + 0, 3 * c + 1),
            &m(3 * r + 0, 3 * c + 2), &m(3 * r + 1, 3 * c + 0),
            &m(3 * r + 1, 3 * c + 1), &m(3 * r + 1, 3 * c + 2),
            &m(3 * r + 2, 3 * c + 0), &m(3 * r + 2, 3 * c + 1),
            &m(3 * r + 2, 3 * c + 2));
      };
 
      for (int gg = 0; gg != nb_gauss_pts; gg++) {
 
        double val = getGaussPts()(2, gg) * 0.5;
 
        auxFun.calculateGriffithForce(1, 1, row_data.getDiffN(gg));
 
        auto t_griffith = get_ftensor_from_vec(auxFun.griffithForce);
        auto t_diff_rho = getFTensor1FromMat<3>(commonData.diffRho);
        for (int rrr = 0; rrr != nb_base_row; rrr++) {
 
          const double a = val * gc * pow(rho, betaGC - 1.) * betaGC;
          FTensor::Tensor2<double, 3, 3> k;
          k(i, j) = t_diff_rho(j) * t_griffith(i) * a;
          // there is no dependence on the other nodes so row == col
          auto tt_mat = get_tensor2(mat, rrr, rrr);
          tt_mat(i, j) += k(i, j);
 
          ++t_diff_rho;
          ++t_griffith;
          ++rho;
        }
      }
 
      int col_nb_dofs = col_data.getIndices().size();
 
      CHKERR MatSetValues(getFEMethod()->snes_B, row_nb_dofs,
                          &*rowIndices.data().begin(), col_nb_dofs,
                          &*col_data.getIndices().data().begin(),
                          &*mat.data().begin(), ADD_VALUES);
 
      MoFEMFunctionReturn(0);
    }
  };
 
  struct MyTriangleFEConstrainsDelta : public MyTriangleFE {
 
    const std::string lambdaFieldName;
    bool ownDeltaVec;
    SmartPetscObj<Vec> deltaVec;
 
    MyTriangleFEConstrainsDelta(MoFEM::Interface &m_field,
                                const std::string &lambda_field_name)
        : MyTriangleFE(m_field), lambdaFieldName(lambda_field_name) {}
 
    MoFEMErrorCode preProcess() {
      MoFEMFunctionBegin;
      Vec l;
      CHKERR VecGhostGetLocalForm(deltaVec, &l);
      double *a;
      CHKERR VecGetArray(l, &a);
      int n;
      CHKERR VecGetLocalSize(l, &n);
      std::fill(&a[0], &a[n], 0);
      CHKERR VecRestoreArray(l, &a);
      CHKERR VecGhostRestoreLocalForm(deltaVec, &l);
      MoFEMFunctionReturn(0);
    }
 
    MoFEMErrorCode postProcess() {
      MoFEMFunctionBegin;
      MOFEM_LOG_CHANNEL("WORLD");
      MOFEM_LOG_FUNCTION();
      MOFEM_LOG_TAG("WORLD", "Arc");
 
      CHKERR VecAssemblyBegin(deltaVec);
      CHKERR VecAssemblyEnd(deltaVec);
      CHKERR VecGhostUpdateBegin(deltaVec, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecGhostUpdateEnd(deltaVec, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecGhostUpdateBegin(deltaVec, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(deltaVec, INSERT_VALUES, SCATTER_FORWARD);
      double nrm;
      CHKERR VecNorm(deltaVec, NORM_2, &nrm);
      MOFEM_LOG_C("WORLD", Sev::noisy, "\tdelta vec nrm  = %9.8e", nrm);
      MoFEMFunctionReturn(0);
    }
  };
 
  struct OpConstrainsDelta
      : public FaceElementForcesAndSourcesCore::UserDataOperator,
        AuxOp {
 
    MoFEM::Interface &mField;
    const std::string lambdaFieldName;
    SmartPetscObj<Vec> &deltaVec;
 
    OpConstrainsDelta(MoFEM::Interface &m_field, int tag, BlockData &block_data,
                      CommonData &common_data,
                      const std::string &lambda_field_name,
                      SmartPetscObj<Vec> &delta_vec)
        : FaceElementForcesAndSourcesCore::UserDataOperator(
              "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
          AuxOp(tag, block_data, common_data), mField(m_field),
          lambdaFieldName(lambda_field_name), deltaVec(delta_vec) {}
 
    AuxFunctions<double> auxFun;
 
    MoFEMErrorCode doWork(int side, EntityType type,
                          DataForcesAndSourcesCore::EntData &data) {
      MoFEMFunctionBegin;
      if (type != MBVERTEX) {
        MoFEMFunctionReturnHot(0);
      }
      // get indices of material DOFs
      CHKERR setIndices(data);
      // get indices of Lagrange multiplier DOFs
      CHKERR setLambdaIndices(this, lambdaFieldName);
      // do calculations
      auxFun.currentCoords.resize(9, false);
      for (int dd = 0; dd != 9; ++dd)
        auxFun.currentCoords[dd] = data.getFieldData()[dd];
 
      CHKERR auxFun.calculateGriffithForce(1, getGaussPts()(2, 0) * 0.5,
                                           data.getDiffN(0));
 
      // assemble area change
      double delta[] = {0, 0, 0};
      for (int nn = 0; nn != 3; ++nn) {
        if (rowLambdaIndices[nn] == -1) {
          if (rowIndices[3 * nn + 0] != -1 || rowIndices[3 * nn + 1] != -1 ||
              rowIndices[3 * nn + 2] != -1) {
            SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                    "Element has no nodes at crack front");
          }
          continue;
        }
        for (int dd = 0; dd != 3; ++dd) {
          int idx = 3 * nn + dd;
          delta[nn] += auxFun.griffithForce[idx] *
                       (data.getFieldData()[idx] - getCoords()[idx]);
        }
      }
      // add nodes to vector of lambda indices
      CHKERR VecSetValues(deltaVec, 3, &rowLambdaIndices[0], delta, ADD_VALUES);
 
      MoFEMFunctionReturn(0);
    }
  };
 
  struct MyTriangleFEConstrains : public MyTriangleFE {
 
    const std::string lambdaFieldName;
    BlockData &blockData;
    SmartPetscObj<Vec> deltaVec;
 
    MyTriangleFEConstrains(MoFEM::Interface &m_field,
                           const std::string &lambda_field_name,
                           BlockData &block_data, SmartPetscObj<Vec> &delta_vec)
        : MyTriangleFE(m_field), lambdaFieldName(lambda_field_name),
          blockData(block_data), deltaVec(delta_vec) {}
 
    MoFEMErrorCode preProcess() {
      MoFEMFunctionBegin;
 
      MOFEM_LOG_CHANNEL("WORLD");
      MOFEM_LOG_FUNCTION();
      MOFEM_LOG_TAG("WORLD", "GriffithConstrain");
 
      const double *delta;
      CHKERR VecGetArrayRead(deltaVec, &delta);
      double sum_of_delta = 0;
      double sum_of_lambda = 0;
      const auto bit_field_number =
          mField.get_field_bit_number(lambdaFieldName);
      switch (snes_ctx) {
      case CTX_SNESSETFUNCTION: {
        for (_IT_NUMEREDDOF_ROW_BY_BITNUMBER_FOR_LOOP_(problemPtr,
                                                       bit_field_number, dit)) {
          if (static_cast<int>(dit->get()->getPart()) ==
              mField.get_comm_rank()) {
            int local_idx = dit->get()->getPetscLocalDofIdx();
            double lambda = dit->get()->getFieldData();
            double W = lambda - blockData.E * delta[local_idx];
            double val = lambda - calMax(W, 0, blockData.r);
            int global_idx = dit->get()->getPetscGlobalDofIdx();
            CHKERR VecSetValue(snes_f, global_idx, val, ADD_VALUES);
          }
        }
      } break;
      case CTX_SNESSETJACOBIAN: {
        for (_IT_NUMEREDDOF_ROW_BY_BITNUMBER_FOR_LOOP_(problemPtr,
                                                       bit_field_number, dit)) {
          if (static_cast<int>(dit->get()->getPart()) ==
              mField.get_comm_rank()) {
            int local_idx = dit->get()->getPetscLocalDofIdx();
            double lambda = dit->get()->getFieldData();
            double W = lambda - blockData.E * delta[local_idx];
            double diffW = 1;
            double val = 1 - diffCalMax_a(W, 0, blockData.r) * diffW;
            int global_idx = dit->get()->getPetscGlobalDofIdx();
            MOFEM_LOG_C("WORLD", Sev::verbose,
                        "Constrains on node %lu diag = %+3.5e "
                        "delta = %+3.5e lambda = %+3.5e",
                        dit->get()->getEnt(), val, delta[local_idx], lambda);
            CHKERR MatSetValue(snes_B, global_idx, global_idx, val, ADD_VALUES);
            sum_of_delta += delta[local_idx];
            sum_of_lambda += lambda;
          }
        }
        MOFEM_LOG_C("WORLD", Sev::noisy,
                    "Sum of delta = %+6.4e Sum of lambda = %+6.4e",
                    sum_of_delta, sum_of_lambda);
 
      } break;
      default:
        break;
      }
      CHKERR VecRestoreArrayRead(deltaVec, &delta);
      MoFEMFunctionReturn(0);
    }
 
    MoFEMErrorCode postProcess() {
      MoFEMFunctionBeginHot;
      MoFEMFunctionReturnHot(0);
    }
  };
 
  struct OpConstrainsRhs
      : public FaceElementForcesAndSourcesCore::UserDataOperator,
        AuxOp {
 
    MoFEM::Interface &mField;
    const std::string lambdaFieldName;
 
    OpConstrainsRhs(MoFEM::Interface &m_field, int tag, BlockData &block_data,
                    CommonData &common_data,
                    const std::string &lambda_field_name)
        : FaceElementForcesAndSourcesCore::UserDataOperator(
              "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
          AuxOp(tag, block_data, common_data), mField(m_field),
          lambdaFieldName(lambda_field_name) {}
 
    AuxFunctions<double> auxFun;
    VectorDouble9 nF;
    // VectorDouble3 nG;
 
    MoFEMErrorCode doWork(int side, EntityType type,
                          DataForcesAndSourcesCore::EntData &data) {
      MoFEMFunctionBegin;
      if (type != MBVERTEX)
        MoFEMFunctionReturnHot(0);
 
      // get indices of material DOFs
      CHKERR setIndices(data);
      CHKERR setLambdaNodes(this, lambdaFieldName);
      CHKERR setLambdaIndices(this, lambdaFieldName);
 
      auxFun.currentCoords.resize(9, false);
      noalias(auxFun.currentCoords) = data.getFieldData();
 
      CHKERR auxFun.calculateGriffithForce(1, getGaussPts()(2, 0) * 0.5,
                                           data.getDiffN(0));
 
      // assemble local vector
      nF.resize(9, false);
      nF.clear();
      for (int nn = 0; nn != 3; nn++) {
        if (rowLambdaIndices[nn] != -1) {
          const auto lambda = lambdaAtNodes[nn];
          for (int dd = 0; dd != 3; dd++) {
            int idx = 3 * nn + dd;
            nF[idx] -= lambda * auxFun.griffithForce[idx];
          }
        }
      }
 
      // assemble the right hand vectors
      CHKERR VecSetValues(getFEMethod()->snes_f, 9, &*rowIndices.data().begin(),
                          &nF[0], ADD_VALUES);
 
      MoFEMFunctionReturn(0);
    }
  };
 
  struct OpConstrainsLhs
      : public FaceElementForcesAndSourcesCore::UserDataOperator,
        AuxOp {
 
    MoFEM::Interface &mField;
    SmartPetscObj<Vec> &deltaVec;
 
    OpConstrainsLhs(MoFEM::Interface &m_field, int tag, BlockData &block_data,
                    CommonData &common_data,
                    const std::string &lambda_field_name,
                    SmartPetscObj<Vec> &delta_vec)
        : FaceElementForcesAndSourcesCore::UserDataOperator(
              "MESH_NODE_POSITIONS", lambda_field_name,
              UserDataOperator::OPROWCOL,
              false // not symmetric operator
              ),
          AuxOp(tag, block_data, common_data), mField(m_field),
          deltaVec(delta_vec) {}
 
    MatrixDouble nG_dX;
    MatrixDouble nF_dLambda;
    MatrixDouble nF_dX;
 
    MatrixDouble jacDelta;
    ublas::vector<double *> jacDeltaPtr;
    VectorDouble d_dA;
    VectorDouble d_Delta;
 
    ublas::vector<adouble> a_Delta;
    AuxFunctions<adouble> a_auxFun;
    ublas::vector<adouble> a_dA;
    ublas::vector<adouble> a_lambdaAtNodes;
 
    MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                          EntityType col_type,
                          DataForcesAndSourcesCore::EntData &row_data,
                          DataForcesAndSourcesCore::EntData &col_data) {
      MoFEMFunctionBegin;
 
      if (row_type != MBVERTEX || col_type != MBVERTEX)
        MoFEMFunctionReturnHot(0);
 
      CHKERR setIndices(row_data);
      CHKERR setLambdaNodes(this, colFieldName);
      CHKERR setLambdaIndices(this, colFieldName);
 
      activeVariables.resize(18, false);
      a_auxFun.referenceCoords.resize(9, false);
      a_auxFun.currentCoords.resize(9, false);
 
      a_Delta.resize(3, false);
      d_Delta.resize(3, false);
      trace_on(tAg);
      for (int dd = 0; dd != 9; dd++) {
        double val = getCoords()[dd];
        a_auxFun.referenceCoords[dd] <<= val;
        activeVariables[dd] = val;
      }
      for (int dd = 0; dd != 9; dd++) {
        double val = row_data.getFieldData()[dd];
        a_auxFun.currentCoords[dd] <<= val;
        activeVariables[9 + dd] = val;
      }
 
      CHKERR a_auxFun.calculateGriffithForce(1, getGaussPts()(2, 0) * 0.5,
                                             row_data.getDiffN(0));
 
      // calculate area change
      std::fill(a_Delta.begin(), a_Delta.end(), 0);
      for (int nn = 0; nn != 3; nn++) {
        if (col_data.getIndices()[nn] == -1)
          continue;
        if (col_data.getIndices()[nn] != rowLambdaIndices[nn]) {
          SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                  "data inconsistency");
        }
        for (int dd = 0; dd != 3; dd++) {
          int idx = 3 * nn + dd;
          if (row_data.getFieldDofs()[idx]->getEnt() !=
              col_data.getFieldDofs()[nn]->getEnt()) {
            SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                    "data inconsistency");
          }
          a_Delta[nn] +=
              a_auxFun.griffithForce[idx] *
              (a_auxFun.currentCoords[idx] - a_auxFun.referenceCoords[idx]);
        }
      }
      for (int nn = 0; nn != 3; nn++) {
        a_Delta[nn] >>= d_Delta[nn];
      }
      trace_off();
 
      jacDelta.resize(3, 18, false);
      jacDeltaPtr.resize(3, false);
      for (int nn = 0; nn != 3; nn++) {
        jacDeltaPtr[nn] = &jacDelta(nn, 0);
      }
 
      {
        int r;
        // play recorder for jacobian
        r = ::jacobian(tAg, 3, 18, &activeVariables[0], &jacDeltaPtr[0]);
        if (r < 3) { // function is locally analytic
          SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
                   "ADOL-C function evaluation with error r = %d", r);
        }
      }
 
      const double *delta;
      Vec delta_vec_local;
      CHKERR VecGhostGetLocalForm(deltaVec, &delta_vec_local);
      int vec_size;
      CHKERR VecGetSize(delta_vec_local, &vec_size);
      CHKERR VecGetArrayRead(delta_vec_local, &delta);
      nG_dX.resize(3, 9, false);
      nG_dX.clear();
      for (int nn = 0; nn != 3; ++nn) {
        int local_idx = col_data.getLocalIndices()[nn];
        if (local_idx != rowLambdaIndicesLocal[nn]) {
          SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                  "data inconsistency");
        }
        if (local_idx == -1)
          continue;
        if (vec_size < local_idx) {
          int g_vec_size;
          CHKERR VecGetLocalSize(deltaVec, &g_vec_size);
          SETERRQ3(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                   "data inconsistency %d < %d (%d)", vec_size, local_idx,
                   g_vec_size);
        }
        double W = col_data.getFieldData()[nn] - blockData.E * delta[local_idx];
        double diff_constrain = -diffCalMax_a(W, 0, blockData.r);
        for (int dd = 0; dd != 9; dd++) {
          double diffW = -blockData.E * jacDelta(nn, 9 + dd);
          nG_dX(nn, dd) = diff_constrain * diffW;
        }
      }
      CHKERR VecRestoreArrayRead(delta_vec_local, &delta);
      CHKERR VecGhostRestoreLocalForm(deltaVec, &delta_vec_local);
      CHKERR MatSetValues(getFEMethod()->snes_B, 3,
                          &*col_data.getIndices().data().begin(), 9,
                          &*row_data.getIndices().data().begin(),
                          &*nG_dX.data().begin(), ADD_VALUES);
 
      // Dervatives of nF
 
      // Linearisation with adol-C
      activeVariables.resize(21, false);
      d_dA.resize(9, false);
      a_dA.resize(9, false);
      a_lambdaAtNodes.resize(3, false);
      a_auxFun.referenceCoords.resize(9, false);
 
      trace_on(tAg);
      for (int dd = 0; dd != 9; dd++) {
        double val = getCoords()[dd];
        a_auxFun.referenceCoords[dd] <<= val;
        activeVariables[dd] = val;
      }
      for (int dd = 0; dd != 9; dd++) {
        double val = row_data.getFieldData()[dd];
        a_auxFun.currentCoords[dd] <<= val;
        activeVariables[9 + dd] = val;
      }
      for (int nn = 0; nn != 3; nn++) {
        double val = col_data.getFieldData()[nn];
        a_lambdaAtNodes[nn] <<= val;
        activeVariables[18 + nn] = val;
      }
 
      CHKERR a_auxFun.calculateGriffithForce(1, getGaussPts()(2, 0) * 0.5,
                                             row_data.getDiffN(0));
 
      std::fill(a_dA.begin(), a_dA.end(), 0);
      for (int nn = 0; nn != 3; nn++) {
        if (col_data.getIndices()[nn] == -1)
          continue;
        for (int dd = 0; dd != 3; dd++) {
          int idx = 3 * nn + dd;
          a_dA[idx] -= a_lambdaAtNodes[nn] * a_auxFun.griffithForce[idx];
        }
      }
      for (int dd = 0; dd != 9; dd++) {
        a_dA[dd] >>= d_dA[dd];
      }
      trace_off();
 
      commonData.tangentGriffithForce.resize(9, 21, false);
      commonData.tangentGriffithForceRowPtr.resize(9);
      for (int dd = 0; dd != 9; ++dd) {
        commonData.tangentGriffithForceRowPtr[dd] =
            &commonData.tangentGriffithForce(dd, 0);
      }
 
      {
        int r;
        // play recorder for jacobian
        r = ::jacobian(tAg, 9, 21, &activeVariables[0],
                       &commonData.tangentGriffithForceRowPtr[0]);
        if (r < 3) { // function is locally analytic
          SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
                   "ADOL-C function evaluation with error r = %d", r);
        }
      }
 
      nF_dLambda.resize(9, 3);
      for (int rr = 0; rr != 9; rr++) {
        for (int cc = 0; cc != 3; cc++) {
          nF_dLambda(rr, cc) = commonData.tangentGriffithForce(rr, 18 + cc);
        }
      }
      CHKERR MatSetValues(getFEMethod()->snes_B, 9, &*rowIndices.data().begin(),
                          3, &*col_data.getIndices().data().begin(),
                          &*nF_dLambda.data().begin(), ADD_VALUES);
 
      nF_dX.resize(9, 9, false);
      for (int rr = 0; rr != 9; rr++) {
        for (int cc = 0; cc != 9; cc++) {
          nF_dX(rr, cc) = commonData.tangentGriffithForce(rr, 9 + cc);
        }
      }
      CHKERR MatSetValues(getFEMethod()->snes_B, 9, &*rowIndices.data().begin(),
                          9, &*row_data.getIndices().data().begin(),
                          &*nF_dX.data().begin(), ADD_VALUES);
 
      MoFEMFunctionReturn(0);
    }
  };
 
  /**
   * \brief arc-length element
   *
   * Calualte residual for arc length method amd vector db, which is derivative
   * of residual.
   *
   */
  struct FrontArcLengthControl : public FEMethod {
 
    int tAg;
    boost::shared_ptr<ArcLengthCtx> arcPtr;
    BlockData &blockData;
    CommonData &commonData;
    std::string lambdaFieldName;
 
    Vec ghostIntLambda;
    double intLambdaArray[2];
    double &intLambdaDelta;
    double &intLambdaLambda;
    double intLambda;
 
    boost::shared_ptr<MyTriangleFE> feIntLambda;
 
    /**
     * \brief assemble internal residual
     * @param block_data        griffith block data
     * @param common_data       griffith common data
     * @param lambda_field_name lagrange multipliers controling crack surface
     * area
     * @param arc_front         reference to FrontArcLengthControl object
     */
    struct OpIntLambda : FaceElementForcesAndSourcesCore::UserDataOperator,
                         GriffithForceElement::AuxOp {
 
      MoFEM::Interface &mField;
      std::string lambdaFieldName;
      FrontArcLengthControl &arcFront;
 
      OpIntLambda(MoFEM::Interface &m_field,
                  GriffithForceElement::BlockData &block_data,
                  GriffithForceElement::CommonData &common_data,
                  const std::string &lambda_field_name,
                  FrontArcLengthControl &arc_front)
          : FaceElementForcesAndSourcesCore::UserDataOperator(
                "MESH_NODE_POSITIONS", lambda_field_name,
                UserDataOperator::OPROW,
                false // not symmetric operator
                ),
            AuxOp(0, block_data, common_data), mField(m_field),
            lambdaFieldName(lambda_field_name), arcFront(arc_front) {}
 
      AuxFunctions<double> auxFun;
 
      MoFEMErrorCode doWork(int side, EntityType type,
                            DataForcesAndSourcesCore::EntData &data) {
        MoFEMFunctionBegin;
        if (type != MBVERTEX) {
          MoFEMFunctionReturnHot(0);
        }
        // get indices of Lagrange multiplier DOFs
        CHKERR setLambdaIndices(this, lambdaFieldName);
        // get dofs increment
        const double *dx;
        CHKERR VecGetArrayRead(arcFront.arcPtr->dx, &dx);
        // get dofs at begining of load step
        const double *x0;
        CHKERR VecGetArrayRead(arcFront.arcPtr->x0, &x0);
        auxFun.referenceCoords.resize(9, false);
        auxFun.currentCoords.resize(9, false);
        for (int dd = 0; dd != 9; dd++) {
          int loc_idx = data.getLocalIndices()[dd];
          if (loc_idx != -1) {
            auxFun.referenceCoords[dd] = x0[loc_idx];
            auxFun.currentCoords[dd] = auxFun.referenceCoords[dd] + dx[loc_idx];
          } else {
            SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Note expected");
          }
        }
        CHKERR VecRestoreArrayRead(arcFront.arcPtr->dx, &dx);
        CHKERR VecRestoreArrayRead(arcFront.arcPtr->x0, &x0);
 
        // calculate crack area increment
        CHKERR auxFun.calculateGriffithForce(1, getGaussPts()(2, 0) * 0.5,
                                             data.getDiffN(0));
 
        // calculate area change
        double delta = 0;
        for (int nn = 0; nn != 3; nn++) {
          if (rowLambdaIndices[nn] == -1)
            continue;
          auto griffith_force =
              getVectorAdaptor(&(auxFun.griffithForce[3 * nn]), 3);
          // griffith_force /= norm_2(griffith_force);
          for (int dd = 0; dd != 3; dd++) {
            const int idx = 3 * nn + dd;
            delta += griffith_force[dd] *
                     (auxFun.currentCoords[idx] - auxFun.referenceCoords[idx]);
          }
        }
        if (delta != delta) {
          SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "nan value");
        }
        arcFront.intLambdaDelta += delta;
        MoFEMFunctionReturn(0);
      };
    };
 
    /**
     * \brief assemble internal residual derivative
     * @param block_data        griffith block data
     * @param common_data       griffith common data
     * @param lambda_field_name lagrange multipliers controling crack surface
     * area
     * @param arc_front         reference to FrontArcLengthControl object
     */
    struct OpAssemble_db : public OpIntLambda {
      int tAg;
      OpAssemble_db(int tag, MoFEM::Interface &m_field,
                    GriffithForceElement::BlockData &block_data,
                    GriffithForceElement::CommonData &common_data,
                    const std::string &lambda_field_name,
                    FrontArcLengthControl &arc_front)
          : OpIntLambda(m_field, block_data, common_data, lambda_field_name,
                        arc_front),
            tAg(tag) {}
 
      AuxFunctions<adouble> a_auxFun;
      VectorDouble gradDelta;
      VectorDouble activeVariables;
 
      MoFEMErrorCode doWork(int side, EntityType type,
                            DataForcesAndSourcesCore::EntData &data) {
        MoFEMFunctionBeginHot;
        if (type != MBVERTEX) {
          MoFEMFunctionReturnHot(0);
        }
        CHKERR setIndices(data);
        // get indices of Lagrange multiplier DOFs
        CHKERR setLambdaIndices(this, lambdaFieldName);
        double d_delta;
        adouble a_delta;
        // get dofs increment
        const double *dx;
        CHKERR VecGetArrayRead(arcFront.arcPtr->dx, &dx);
        // get dofs at begining of load step
        const double *x0;
        CHKERR VecGetArrayRead(arcFront.arcPtr->x0, &x0);
        a_auxFun.referenceCoords.resize(9, false);
        a_auxFun.currentCoords.resize(9, false);
 
        activeVariables.resize(18, false);
        trace_on(tAg);
        for (int dd = 0; dd != 9; dd++) {
          int loc_idx = data.getLocalIndices()[dd];
          if (loc_idx != -1) {
            a_auxFun.referenceCoords[dd] <<= x0[loc_idx];
            activeVariables[dd] = x0[loc_idx];
          } else {
            SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Note expected");
          }
        }
        for (int dd = 0; dd != 9; dd++) {
          int loc_idx = data.getLocalIndices()[dd];
          if (loc_idx != -1) {
            a_auxFun.currentCoords[dd] <<= x0[dd] + dx[loc_idx];
            activeVariables[9 + dd] = x0[loc_idx] + dx[loc_idx];
          } else {
            SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Note expected");
          }
        }
        // calculate crack area increment
 
        CHKERR a_auxFun.calculateGriffithForce(1, getGaussPts()(2, 0) * 0.5,
                                               data.getDiffN(0));
        // calculate area change
        a_delta = 0;
        for (int nn = 0; nn != 3; nn++) {
          if (rowLambdaIndices[nn] == -1)
            continue;
          auto griffith_force =
              getVectorAdaptor(&(a_auxFun.griffithForce[3 * nn]), 3);
          // griffith_force /= sqrt(griffith_force[0] * griffith_force[0] +
          //                         griffith_force[1] * griffith_force[1] +
          //                         griffith_force[2] * griffith_force[2]);
          for (int dd = 0; dd != 3; dd++) {
            const int idx = 3 * nn + dd;
            a_delta += griffith_force[dd] * (a_auxFun.currentCoords[idx] -
                                             a_auxFun.referenceCoords[idx]);
          }
        }
        a_delta >>= d_delta;
        trace_off();
        CHKERR VecRestoreArrayRead(arcFront.arcPtr->dx, &dx);
        CHKERR VecRestoreArrayRead(arcFront.arcPtr->x0, &x0);
 
        gradDelta.resize(18, false);
        {
          int r;
          // play recorder for jacobian
          r = ::gradient(tAg, 18, &activeVariables[0], &gradDelta[0]);
          if (r < 3) { // function is locally analytic
            SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
                     "ADOL-C function evaluation with error r = %d", r);
          }
        }
 
        for (int dd = 0; dd != 9; dd++) {
          double val = arcFront.arcPtr->alpha * gradDelta[9 + dd];
          CHKERR VecSetValue(arcFront.arcPtr->db, data.getIndices()[dd], val,
                             ADD_VALUES);
        }
 
        MoFEMFunctionReturnHot(0);
      }
    };
 
    FrontArcLengthControl(int tag, GriffithForceElement::BlockData &block_data,
                          GriffithForceElement::CommonData &common_data,
                          const std::string &lambda_field_name,
                          boost::shared_ptr<ArcLengthCtx> &arc_ptr)
        : tAg(tag), arcPtr(arc_ptr), blockData(block_data),
          commonData(common_data), lambdaFieldName(lambda_field_name),
          intLambdaDelta(intLambdaArray[0]),
          intLambdaLambda(intLambdaArray[1]) {
      int two[] = {0, 1};
      int nb_local = arcPtr->mField.get_comm_rank() == 0 ? 2 : 0;
      int nb_ghost = nb_local ? 0 : 2;
      ierr = VecCreateGhostWithArray(PETSC_COMM_WORLD, nb_local, 2, nb_ghost,
                                     two, intLambdaArray, &ghostIntLambda);
      CHKERRABORT(PETSC_COMM_WORLD, ierr);
      feIntLambda = boost::shared_ptr<MyTriangleFE>(
          new GriffithForceElement::MyTriangleFE(arcPtr->mField));
    }
    virtual ~FrontArcLengthControl() {
      ierr = VecDestroy(&ghostIntLambda);
      CHKERRABORT(PETSC_COMM_WORLD, ierr);
    }
 
    MoFEMErrorCode preProcess() {
      MoFEMFunctionBegin;
      switch (snes_ctx) {
      case CTX_SNESSETFUNCTION: {
        CHKERR VecAssemblyBegin(snes_f);
        CHKERR VecAssemblyEnd(snes_f);
        CHKERR calculateDxAndDlambda(snes_x);
        CHKERR calculateDb();
      } break;
      default:
        break;
      }
      MoFEMFunctionReturn(0);
    }
 
    /** \brief Calculate internal lambda
     */
    double calculateLambdaInt() {
      MoFEMFunctionBeginHot;
      return intLambda;
      MoFEMFunctionReturnHot(0);
    }
 
    MoFEMErrorCode operator()() {
      MoFEMFunctionBegin;
      MOFEM_LOG_CHANNEL("WORLD");
      MOFEM_LOG_FUNCTION();
      MOFEM_LOG_TAG("WORLD", "Arc");
 
      switch (snes_ctx) {
      case CTX_SNESSETFUNCTION: {
        arcPtr->res_lambda = calculateLambdaInt() - arcPtr->alpha * arcPtr->s;
        CHKERR VecSetValue(snes_f, arcPtr->getPetscGlobalDofIdx(),
                           arcPtr->res_lambda, ADD_VALUES);
      } break;
      case CTX_SNESSETJACOBIAN: {
        arcPtr->dIag = 2 * arcPtr->dLambda * arcPtr->alpha *
                       pow(arcPtr->beta, 2) * arcPtr->F_lambda2;
        MOFEM_LOG_C("WORLD", Sev::noisy, "diag = %9.8e", arcPtr->dIag);
        CHKERR MatSetValue(snes_B, arcPtr->getPetscGlobalDofIdx(),
                           arcPtr->getPetscGlobalDofIdx(), 1, ADD_VALUES);
      } break;
      default:
        break;
      }
      MoFEMFunctionReturn(0);
    }
 
    MoFEMErrorCode postProcess() {
      MoFEMFunctionBegin;
      switch (snes_ctx) {
      case CTX_SNESSETFUNCTION: {
        CHKERR VecAssemblyBegin(snes_f);
        CHKERR VecAssemblyEnd(snes_f);
      } break;
      case CTX_SNESSETJACOBIAN: {
        CHKERR VecGhostUpdateBegin(arcPtr->ghostDiag, INSERT_VALUES,
                                   SCATTER_FORWARD);
        CHKERR VecGhostUpdateEnd(arcPtr->ghostDiag, INSERT_VALUES,
                                 SCATTER_FORWARD);
      } break;
      default:
        break;
      }
      MoFEMFunctionReturn(0);
    }
 
    /** \brief Calculate db
     */
    MoFEMErrorCode calculateDb() {
      MoFEMFunctionBegin;
 
      MOFEM_LOG_CHANNEL("WORLD");
      MOFEM_LOG_FUNCTION();
      MOFEM_LOG_TAG("WORLD", "Arc");
 
      double alpha = arcPtr->alpha;
      double beta = arcPtr->beta;
      double d_lambda = arcPtr->dLambda;
 
      CHKERR VecZeroEntries(ghostIntLambda);
      CHKERR VecGhostUpdateBegin(ghostIntLambda, INSERT_VALUES,
                                 SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(ghostIntLambda, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecZeroEntries(arcPtr->db);
      CHKERR VecGhostUpdateBegin(arcPtr->db, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(arcPtr->db, INSERT_VALUES, SCATTER_FORWARD);
 
      // calculate internal lambda
      feIntLambda->getOpPtrVector().clear();
      feIntLambda->getOpPtrVector().push_back(new OpIntLambda(
          arcPtr->mField, blockData, commonData, lambdaFieldName, *this));
      feIntLambda->getOpPtrVector().push_back(new OpAssemble_db(
          tAg, arcPtr->mField, blockData, commonData, lambdaFieldName, *this));
      CHKERR arcPtr->mField.loop_finite_elements(
          problemPtr->getName(), "CRACKFRONT_AREA_ELEM", (*feIntLambda));
      const double *dx;
      const auto bit_field_number = getFieldBitNumber(lambdaFieldName);
      CHKERR VecGetArrayRead(arcPtr->dx, &dx);
      for (_IT_NUMEREDDOF_ROW_BY_BITNUMBER_FOR_LOOP_(problemPtr,
                                                     bit_field_number, dit)) {
        if (static_cast<int>(dit->get()->getPart()) !=
            arcPtr->mField.get_comm_rank())
          continue;
        int idx = dit->get()->getPetscLocalDofIdx();
        double val = dx[idx];
        CHKERR VecSetValue(ghostIntLambda, 1, val, ADD_VALUES);
      }
      CHKERR VecRestoreArrayRead(arcPtr->dx, &dx);
 
      CHKERR VecAssemblyBegin(ghostIntLambda);
      CHKERR VecAssemblyEnd(ghostIntLambda);
      CHKERR VecGhostUpdateBegin(ghostIntLambda, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecGhostUpdateEnd(ghostIntLambda, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecGhostUpdateBegin(ghostIntLambda, INSERT_VALUES,
                                 SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(ghostIntLambda, INSERT_VALUES, SCATTER_FORWARD);
 
      CHKERR VecAssemblyBegin(arcPtr->db);
      CHKERR VecAssemblyEnd(arcPtr->db);
      CHKERR VecGhostUpdateBegin(arcPtr->db, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecGhostUpdateEnd(arcPtr->db, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecGhostUpdateBegin(arcPtr->db, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(arcPtr->db, INSERT_VALUES, SCATTER_FORWARD);
 
      intLambda = alpha * intLambdaDelta +
                  alpha * beta * beta * d_lambda * d_lambda * arcPtr->F_lambda2;
 
      MOFEM_LOG_C(
          "WORLD", Sev::noisy,
          "\tintLambdaLambda = %9.8e intLambdaDelta = %9.8e intLambda = %9.8e",
          intLambdaLambda, intLambdaDelta, intLambda);
 
      double nrm;
      CHKERR VecNorm(arcPtr->db, NORM_2, &nrm);
      MOFEM_LOG_C("WORLD", Sev::noisy, "\tdb nrm  = %9.8e", nrm);
 
      MoFEMFunctionReturn(0);
    }
 
    MoFEMErrorCode calculateDxAndDlambda(Vec x) {
      MoFEMFunctionBegin;
      MOFEM_LOG_CHANNEL("WORLD");
      MOFEM_LOG_FUNCTION();
      MOFEM_LOG_TAG("WORLD", "Arc");
 
      // Calculate dx
      CHKERR VecGhostUpdateBegin(arcPtr->x0, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(arcPtr->x0, INSERT_VALUES, SCATTER_FORWARD);
 
      Vec l_x, l_x0, l_dx;
      CHKERR VecGhostGetLocalForm(x, &l_x);
      CHKERR VecGhostGetLocalForm(arcPtr->x0, &l_x0);
      CHKERR VecGhostGetLocalForm(arcPtr->dx, &l_dx);
      {
        double *x_array, *x0_array, *dx_array;
        CHKERR VecGetArray(l_x, &x_array);
        CHKERR VecGetArray(l_x0, &x0_array);
        CHKERR VecGetArray(l_dx, &dx_array);
        int size =
            problemPtr->getNbLocalDofsRow() + problemPtr->getNbGhostDofsRow();
        for (int i = 0; i != size; ++i) {
          dx_array[i] = x_array[i] - x0_array[i];
        }
        CHKERR VecRestoreArray(l_x, &x_array);
        CHKERR VecRestoreArray(l_x0, &x0_array);
        CHKERR VecRestoreArray(l_dx, &dx_array);
      }
      CHKERR VecGhostRestoreLocalForm(x, &l_x);
      CHKERR VecGhostRestoreLocalForm(arcPtr->x0, &l_x0);
      CHKERR VecGhostRestoreLocalForm(arcPtr->dx, &l_dx);
 
      // Calculate dlambda
      if (arcPtr->getPetscLocalDofIdx() != -1) {
        double *array;
        CHKERR VecGetArray(arcPtr->dx, &array);
        arcPtr->dLambda = array[arcPtr->getPetscLocalDofIdx()];
        array[arcPtr->getPetscLocalDofIdx()] = 0;
        CHKERR VecRestoreArray(arcPtr->dx, &array);
      }
      CHKERR VecGhostUpdateBegin(arcPtr->ghosTdLambda, INSERT_VALUES,
                                 SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(arcPtr->ghosTdLambda, INSERT_VALUES,
                               SCATTER_FORWARD);
 
      // Calculate dx2
      double x_nrm, x0_nrm;
      CHKERR VecNorm(x, NORM_2, &x_nrm);
      CHKERR VecNorm(arcPtr->x0, NORM_2, &x0_nrm);
      CHKERR VecDot(arcPtr->dx, arcPtr->dx, &arcPtr->dx2);
      MOFEM_LOG_C("WORLD", Sev::noisy,
                  "\tx norm = %6.4e x0 norm = %6.4e dx2 = %6.4e", x_nrm, x0_nrm,
                  arcPtr->dx2);
      MoFEMFunctionReturn(0);
    }
  };
};
 
} // namespace FractureMechanics
 
#endif //__GRIFFITH_FORCE_ELEMENT_HPP__